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THE  SECOND  ZIENKIEWICZ  LECTURE 


ERROR  ESTIMATION  AND  CONTROL 
IN 

COMPUTATIONAL  FLUID  DYNAMICS 
by 

J  Tinsley  Oden 


This  paper  describes  theory  and  methods  for  developing  a  posteriori  error  estimates 
and  an  adaptive  strategy  for  /ip-finite  element  approximations  of  the  incompressible 
Navier-Stokes  equations. 

For  an  error  estimation,  use  is  made  of  a  new  approach  which  is  based  on  the  work 
of  Ainsworth,  Wu  and  the  author.  That  theory  has  been  shown  to  produce  good 
results  for  general  elliptic  systems  and  general  ftp-finite  element  methods.  Recently, 
these  techniques  have  been  extended  to  the  Navier-Stokes  equations 

A  three-step  adaptive  algorithm  is  also  described  which  produces  reasonable  hp 
meshes  very  efficiently.  These  techniques  are  applied  to  representative  transient  and 
steady  state  problems  of  incompressible  viscous  flow. 
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Reliability  of  computational  mechanics 

The  aim  of  computational  mechanics  is  to  provide  a  reliable  description  of  physical 
phenomena  related  to  the  aim  of  the  analysis.  Typically,  for  example,  for  preventing  a 
failure,  avoiding  various  engineering  experiments,  etc. 

The  framework  of  computational  analysis  is  a)  mathematical  modelling;  b)  numerical 
solution;  and  c)  interpretation  of  the  results.  The  errors  of  all  these  parts  have  to  be 
under  control  to  achieve  the  desired  reliability. 

The  talk  will  address 

a)  The  modelling  of  plasticity  problems  under  cyclic  loadiing  for  5454  aluminium 
alloy  in  the  H-34  temper. 

b)  The  adaptive  modelling  of  the  heat  problem  in  thin  domains  and  a-posteriori 
estimation  of  the  modelling  errors. 

c)  The  task  of  a-posteriori  error  estimation  in  finite  elements  and  the  theoretical 
principles  of  comparing  today's  major  estimators  used  in  the  literature;  and  the 
concrete  results  related  to  the  performance  of  the  estimators  in  the  case  of  nearly 
patchwise  uniform  meshes. 


H  BERGER,  G  WARNECKE  and  W  L  WENDLAND 
Coupling  of  FEM  and  BEM  for  transonic  flows 

Consider  a  two-dimensional  stationary  compressible  inviscid  transonic  flow  around  a 
given  airfoil.  As  a  mathematical  model  we  use  the  full  potential  equation  of 
isoenergetic,  homentropic  flows  of  an  ideal  gas.  As  is  well  known,  this  equation 
degenerates  where  the  local  density  vanishes  and  for  subsonic  travelling  velocities  the 
flow  develops  regions  with  supersonic  velocity  separated  by  down-stream  shocks  from  the 
down-stream  subsonic  flow  region.  For  the  numerical  computation  we  choose  a  bounded 
ring  domain  around  the  airfoil  and  use  the  generalized  finite  element  least  squares 
method  developed  by  Glowinsky  et  al. 

The  flow  must  contain  an  appropriate  vorticity  and  the  jump  of  the  potential  along  a 
slit  generated  by  the  trailing  edge  subject  to  the  Kutta-Joukowski  condition.  In  order 
to  avoid  nonphysical  solutions  we  incorporate  an  appropriate  penalty  term,  in  addition 
to  the  least  squares  functional,  which  provides  the  right  selection  principle. 
Traditionally,  the  exterior  Neumann  condition  on  the  boundary  of  the  computational 
domain  is  given  by  the  flux  generated  by  an  exterior  parallel  flow.  However,  if  the 
exterior  flow  was  known,  the  flux  would  be  given  by  the  corresponding  cononna! 
derivative  of  the  exterior  potential  associated  with  the  continuous  Diricblet  boundary 
data  belonging  to  the  interior  potential  field.  The  mapping  of  the  Dirichlet  data  into  the 
exterior  flux  is  given  by  the  nonlocal  Steklov-Poincarf  operator  associated  with  the 
exterior  field.  Replacing  the  exterior  flow  by  a  linear  Prandtl-Glavert  equation  given  by 
the  quantities  of  the  parallel  flow  at  infinity,  the  Steklov-Poincari  operator  can  explicitly 
be  defined  in  terms  of  appropriate  boundary  integral  operators  defined  on  the  artificial 


2 


boundary  of  the  bounded  computational  ring  domain  around  the  airfoil.  Tbe  numerical 
approximation  of  the  Steklov-Poincare  operator  then  can  be  given  in  terms  of  boundary 
element  formulations.  If  we  assume  that  the  exterior  model  flow  problem  admits  a 
unique  solution  satisfying  the  Kutta-Joukowski  condition  at  the  trailing  edge  having  there 
a  subsonic  local  flow  region  and  under  the  assumption  that  the  flow  speed  through  the 
artificial  boundary  is  subsonic,  we  can  show  that  a  corresponding  coupled  finite  element  - 
boundary  element  optimization  problem  provides  a  sequence  of  approximate  solutions 
which  converges  to  the  unique  solution  of  an  exterior  transmission  problem  defined  by 
the  weak  solution  of  the  transonic  full  potential  equation  in  the  ring  domain  and  the 
exterior  Prandtl-Glavert  equation.  The  convergence  proof  is  essentially  based  on  a 
compactness  result  by  Mandel  and  Neias  in  connection  with  standard  approximation 
results  for  finite  elements  and  for  boundary  elements.  Tbe  computational  results  are  in 
good  agreement  with  the  aforementioned  analysis. 


G  BUGEDA  and  E  ONATE 

A  study  of  mesh  optimality  criteria  in  adaptive  finite  element  anlysis 

The  evaluation  of  discretization  error  and  the  design  of  suitable  meshes  via  adaptive 
mesh  refinement  (AMR)  are  nowadays  two  of  the  challenging  issues  in  the  finite  element 
method  (FEM). 

in  this  paper  a  methodology  for  deriving  AMR  procedures  is  proposed.  The  basis  of 
the  approach  is  the  decoupling  of  the  concepts  of  error  measure  and  mesh  optimality 
criteria.  This  leads  to  the  definitions  of  global  and  local  error  parameters  from  which 
the  element  refinement  strategy  can  be  simply  ob'ained.  Moreover,  this  allows  us  to 
identify  clearly  the  convergence  rates  of  the  global  and  local  error  norms,  which  strongly 
influence  the  expressions  of  the  element  refinement  parameter.  It  is  shown  that  an 
inaccurate  evaluation  of  this  important  parameter  can  lead  to  oscillations  in  the 
refinement  process. 

The  methodology  proposed  is  particularized  for  two  mesh  optimality  criteria  using 
Zienkiewicz-Zhu  error  estimator  (1),  First  the  standard  criterion  of  equal  distribution 
of  the  global  error  over  all  the  elements  is  studied.  We  will  show  here  that  a  careful 
interpretation  of  the  concepts  of  global  error  and  optimal  mesh  leads  in  this  case  to  an 
enhanced  expression  of  the  element  refinement  parameter  slightly  different  from  that 
typically  used  in  literature. 

The  second  mesb  optimality  criterion  studied  is  based  in  tbe  equal  distribution  of  the 
error  per  unit  area  or  volume  (i.e.  tbe  specific  error).  This  strategy  allows  us  to 
concentrate  more  and  smaller  elements  in  zones  where  the  gradients  of  the  problem 
unknowns  (i.e.  stresses  in  structural  problems)  are  higher,  as  should  be  expected  from 
the  engineering  point  of  view. 

The  layout  of  the  paper  is  the  following.  First  the  concepts  of  error  estimator  and 
acceptable  solution  are  briefly  revisited.  Then  a  general  expression  for  the  element 
refinement  parameter  in  terms  of  global  and  local  error  parameters  is  derived.  The 
mesh  optimality  criteria  based  on  global  and  specific  error  distributions  are  then 
analyzed  and  the  expressions  for  the  local  and  global  error  parameters  and  that  of  tbe 
element  refinement  parameter  are  obtained  for  each  case.  Finally,  the  efficiency  of  the 
two  AMR  methodologies  studied  is  compared  in  tbe  analysis  of  plane  strain,  plate,  shell 
and  potential  flow  problems. 


3 


References: 

Zienkiewicz,  O  C  and  Zhu,  J  Z,  A  simple  error  estimator  and  adaptive  procedure  for 
practical  engineering  analysis.  Int.  Num.  Meth.  Engrg.  24,  337-357,  1987. 


M  J  CROCHET 

Finite  elements  for  polymer  flow:  A  comparison 

Over  the  last  five  years,  several  finite  element  methods  have  been  introduced  for 
analyzing  viscoelastic  flow.  In  the  present  paper,  we  examine  a  class  of  finite  elements 
of  the  mixed  type  where  the  stress  and  velocity  components  and  the  pressure  are 
independent  variables;  they  are  combined  with  various  methods  of  upwinding  for  the 
constitutive  equations.  We  have  selected  several  benchmark  problems;  with  and  without 
geometrical  singularities,  with  and  without  change  of  type  of  the  vorticity  equation.  We 
show  that  it  is  presently  difficult  to  select  "the"  optimal  method  for  all  flows  although 
some  algorithms  perform  much  better  than  others.  We  also  explain  some  outstanding 
problems  for  the  simulation  of  viscoelastic  flow. 


C  JOHNSON 

Adaptive  finite  element  methods  in  computational  mechanics 

We  give  an  overview  of  our  recent  work  on  adaptive  finite  element  methods  for 
problems  in  solid  and  fluid  mechanics  including  elasticity  and  plasticity  problems,  and 
incompressible  and  compressible  fluid  flow.  We  present  new  a-posteriori  error  estimates 
for  the  indicated  classes  of  problems  and  show  that  the  corresponding  adaptive 
algorithms  for  quantitative  error  control  are  reliable  and  efficient.  We  discuss  the 
theoretical  principles  underlying  the  a-posteriori  error  estimates  and  present  the  results 
of  numerical  experiments. 


V  MAZ’YA 

Approximate  approximations 

The  starting  idea  is  to  represent  an  arbitrary  function  as  a  linear  combination  of  basis 
functions,  which,  in  contrast  to  splines,  form  an  approximate  partition  of  unity.  The 
approximations  obtained  do  not  converge  as  the  mesh  size  tends  to  zero.  The  lack  of 
the  convergence  is  compensated  for,  first  of  all,  by  the  flexibility  in  the  choice  of  the 
basis  functions  and  by  the  simplicity  of  the  generalization  to  the  multidimensional  case. 
Another,  and  probably  the  most  important  advantage,  is  the  possibility  to  obtain  explicit 
formulae  for  values  of  various  integral  and  pseudodifferential  operators  of  mathematical 
physics  applied  to  the  basis  functions. 

This  property,  when  used  for  solving  elliptic  boundary  value  problems,  leads  to  a  new 
approach  to  the  discretization  of  boundary  integral  equations,  which  is  not  based  upon 
the  decomposition  of  the  boundary  into  elements.  In  this  approach  the  coefficients  of 
the  resulting  algebraic  system  depend  only  on  the  coordinates  of  a  finite  number  of 
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boundary  points;  hence  the  name  Boundary  Point  Method  (BPM)  seems  quite  natural. 

The  first  step  of  the  BPM  consists  of  the  approximation  of  the  boundary  data  and  of 
the  potential  densities  by  a  linear  combination  of  basis  functions,  each  of  them  being 
either  concentrated  near  a  particular  boundary  point  or  decreasing  rapidly  with  the 
distance  from  this  point.  In  the  second  step  the  calculation  of  the  potentials,  whose 
densities  are  the  basis  functions,  is  carried  out.  I  discuss  only  one  concrete  variant  of 
the  BPM  which  leads  to  an  algebraic  system  with  coefficients  that  are  explicitly  expressed 
-  up  to  a  simple  quadrature. 

The  approximation  method  can  be  employed  also  to  solve  various  non- stationary 
problems.  Its  effectiveness  is  demonstrated  by  the  numerical  solution  of  Cauchy  problem 
for  quasilinear  parabolic  equations.  An  important  feature  of  the  algorithms  is  that  they 
are  both  explicit  and  stable  under  much  milder  restrictions  to  the  time  step,  depending 
on  the  size  of  the  space  grid,  in  comparison  with  usual  explicit  difference  schemes. 


K  W  MORTON  and  E  SLILI 

A  priori  and  a  posteriori  error  analysis  of  finite  volume  methods 

Over  the  last  twenty  years,  finite  volume  methods  have  enjoyed  great  popularity  in  the 
computational  aerodynamics  community  and  have  been  widely  used  for  the  numerical 
simulation  of  transonic  and  supersonic  flows.  The  basic  idea  behind  their  construction 
is  to  exploit  the  divergence  form  of  the  conservation  laws  governing  the  flows  by 
integrating  them  over  finite  volumes,  and  converting  the  volume  integrals  into  contour 
integrals  which  are  then  discretised. 

However,  despite  their  practical  success,  their  theoretical  foundation  is  still 
unsatisfactory  and  their  stability  and  accuracy  properties  are  not  well  undestood. 

The  aiiii  of  this  lecture  is  to  present  an  overview  of  some  recent  theoretical 
developments  concerning  finite  volume  approximations  of  hyperbolic  problems,  and 
parabolic  problems,  with  dominant  hyperbolic  behaviour.  Operating  within  the 
framework  of  non-conforming  Petrov-Galerkin  methods,  we  develop  the  analysis  with 
particular  emphasis  on  the  cell-vertex  scheme. 

A  key  ingredient  of  the  a  priori  error  analysis  is  a  discrete  Girding  inequality  which 
ensures  coercivity  in  a  generalised  sense.  Hence  we  derive  optimal  error  bounds  in 
various  mesh-dependent  norms  and  investigate  the  effects  of  mesh  distortion  on  the 
accuracy  of  the  scheme. 

Through  a  decomposition  of  the  global  error  into  its  local  and  propagating  part,  we 
derive  residual-based  local  a  posteriori  error  estimates  and  indicate  their  relevance  in 
adaptive  mesh-refinement. 

Our  theoretical  results  will  be  illustrated  by  numerical  experiments. 


S  SHAW,  M  K  WARBY  and  J  R  WHITEMAN 

Numerical  techniques  for  problems  of  quasistatic  and  dynamic  viscoelasticity 

The  problem  of  linear  viscoelastic  stress  analysis  of  solid  bodies  is  defined  in  terms  of 
the  governing  momentum  equation  and  two  alternative  hereditary  constitutive  laws. 
Initially  it  is  assumed  that  the  deformation  takes  place  under  quasistatic  conditions  and 
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the  inertial  term  in  the  momentum  equation  is  neglected,  providing  an  equilibrium 
problem  at  each  time  level.  Weak  formulations  of  this  problem  are  given  for  each  of  the 
constitutive  relationships.  These  are  (semi)  discretised  in  time  using,  variously, 
numerical  quadrature  and  finite  difference  replacements  to  produce  the  two  fully 
discrete  formulations,  and  error  estimates  for  the  solutions  of  these  in  the  H'-norm  are 
stated.  Two  independent  numerical  algorithms  result  from  the  fully  discrete 
formulations.  Numerical  results  arising  from  these  are  illustrated  for  model  problems 
and  for  a  problem  involving  a  nylon  material.  These  results  illustrate  the  range  of 
loadings  and  times  over  which  the  linear  viscoelastic  constitutive  models  are  valid.  Some 
discussion  of  extensions  to  nonlinear  models  is  given. 

Finally  the  case  where  the  neglecting  of  the  acceleration  term  in  the  momentum 
equation  is  not  justified,  is  discussed.  In  this  context  the  formulation  is  in  terms  of  an 
integrodifferential  equation  of  either  hyperbolic  or  parabolic  type,  depending  on  which 
constitutive  law  is  adopted.  For  each  of  these  numerical  schemes  are  presented  together 
with  numerical  results. 


Mary  F  WHEELER  and  N  HARDING 

A  characteristics-mixed  method  for  modelling  bioremediation  in  groundwater 

Transport  of  contaminants  in  groundwater  is  frequently  described  by  advective-diffusion 
equations.  Numerical  procedures  based  on  the  coupling  of  characteristic  methods  with 
Galerkin  finite  elements  or  point  centered  finite  difference  methods  for  treating  the 
diffusion  (Eulerian  Lagrangian,  modified  method  of  characteristics  (MMOC).  ELLAM) 
have  been  applied  with  much  success.  One  of  the  main  advantages  of  these  approaches 
is  that  large  time  steps  may  be  employed  for  a  given  accuracy.  Disadvantages  of  some 
of  these  schemes  include  difficulties  of  implementation  in  three  dimensions  and 
nonconserva  lion . 

Here  we  formulate  and  analyze  in  three  spatial  dimensions  a  conservative 
characteristics  mixed  algorithm  for  linear  advection  problems.  In  this  procedure  the  time 
derivative  and  the  advection  terms  are  combined  as  a  directional  derivative  and  the 
spatial  derivatives  are  treated  as  in  a  mixed  finite  element  method.  The  concentration 
approximating  space  is  assumed  to  contain  piecewise  discontinuous  constants  defined 
over  the  triangulation. 

In  our  presentation  we  will  discuss  application  to  a  bioremediation  problem  involving 
a  three  dimensional  six  component  model.  We  have  employed  the  lowest  order  Raviart- 
Thomas  spaces  in  our  numerical  experiments.  Concentrations  and  their  dispersive  fluxes 
are  obtained  and  the  latter  are  used  to  construct  a  higher  order  approximation 
(discontinuous  trilinears)  for  concentration  in  each  cell.  A  slope-limited  scheme  is 
employed  when  necessary  to  prevent  overshoot  and  undershoot.  The  kinetics  are  treated 
by  operator  splitting. 
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ALL  PAPERS  IN  LECTURE  SESSIONS 
J  AALTO,  J  LOUHIVIRTA  and  J  R  WHITEMAN 
On  gradient  smoothing  in  grids  with  singularity  element  patches 

One  possibility  for  dealing  with  corner  singularities  in  the  finite  element  analysis  of 
boundary  value  problems  is  to  use  a  special  subgrid  of  curved  elements,  called  a 
singularity  element  patch,  around  each  corner  point.  This  is  done  by  introducing  a 
suitable  geometrical  mapping  between  the  parametric  u,v-plane  and  the  physical  x,y- 
plane.  In  the  u,v-plane  the  geometry  of  the  elements  is  linear.  In  the  x,y-plane  the 
elements  are  distorted  according  to  the  equations  of  mapping. 

This  paper  shows  how  typical  gradient  smoothing  techniques  can  be  modified  so  that 
they  work  properly  in  grids  with  singularity  element  patches.  The  idea  of  the 
modification  is  based  on  reasonable  smoothed  approximation  for  the  gradient  g  =  V<p 
Within  the  elements  outside  the  singularity  element  patches  the  components  g,  =  dift/dx 
and  gy  =  d<p/i)y  of  the  gradient  are  represented  typically  as 

iM  *  52  •  ij  m  ty**  ■  0) 

>•1  i*i 

In  the  vicinity  of  the  singularity  point,  however,  the  gradient  components  g,  and  g, 
are  singular.  Therefore  it  is  not  reasonable  to  use  approximations  (1),  as  such,  within 
the  elements  inside  a  singularity  element  patch.  Instead  of  these,  the  components  g. 
=  d<f>/d u  and  g,  s  d<f>/dv,  which  are  regular,  are  represented  as 

*.  ■  £  Nt«*  •  (2) 
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In  terms  of  these  the  components  g,  and  gy  are  obtained  using  known  relations  between 
the  gradient  components  g„  gy  and  g,,  g,. 

Such  modifications  for  three  well  known  smoothing  procedures:  Nodal  Averaging 
(N  A),  Global  Least  Squares  smoothing  (GLS)  and  the  Zienkiewicz-Zbu  Superconvergent 
Patch  Recovery  technique  (SPR)  are  dealt  with.  Numerical  examples  of  both  (a)  the 
Poisson’s  equation  and  (b)  the  plane  elasticity  problem  are  given. 


R  T  ACKROYD,  C  R  E  de  OLIVEIRA  and  A  J  H  GODDARD 

Time-dependent  finite  element  solution  of  the  Boltzmann  equation  for  photon  transport 

A  method  of  solving  the  time-dependent  Boltzmann  equation  for  particle  transport  in 
the  case  of  material-radiation  intercation  is  described.  Spatial  dependence  of  the 
solution  is  treated  by  the  finite  element  method.  A  discontinuous  variational  formulation 
is  employed  in  order  to  cope  with  the  abrupt  changes  of  the  physical  characteristics  of 
the  system  with  time.  Particle  scattering  is  treated  by  means  of  a  spherical  harmonic 
expansion  for  particle  direction.  The  method  is  illustrated  by  the  solution  of  the  non- 


7 


linear  equation  for  photon  transport  in  a  reactive  atmosphere  by  presenting  solutions  for 
the  classical  Marshak  thermal  wave  problem.  The  results  show  that  the  method  is 
capable  of  accurately  predicting  the  correct  speed  of  the  thermal  wave  even  in  the  case 
of  very  coarse  meshes. 


R  L  ACT1S  and  B  A  SZAB6 

» 

Stresses  computed  from  hierarchic  plate  models 

Stress  distributions  near  boundaries,  discontinuities  and,  in  the  case  of  laminated  plates, 
at  interfaces  are  generally  very  different  from  the  stress  distribution  in  the  interior 
regions  of  structural  plates.  Boundary  layer  effects  are  normally  present,  and  the 
problem  in  those  regions  is  essentially  three-dimensional.  Hiearchic  models  for 
laminated  plates  make  it  possible  to  approximate  the  three-dimensional  problem  through 
the  solution  of  two-dimensional  problems  without  the  expense  of  a  fully  three- 
dimensional  analysis. 

Higher  order  models  can  be  selected  in  regions  of  local  discontinuities  or  near 
boundaries,  while  low-order  models  are  used  in  the  smooth  interior  regions  of  the 
domain.  Improved  transverse  stresses  can  be  obtained  for  any  model  of  the  hierarchy 
by  using  an  extraction  technique.  In  particular,  the  integration  of  the  equilibrium 
equations  has  shown  substantial  improvement  over  the  direct  compulation  for  the 
transverse  shear  and  normal  stresses.  Examples  are  presented  which  show  that 
hierarchic  models  are  capable  of  approximating  the  fully  three-dimensional  solution. 


M  AINSWORTH 

A  posteriori  error  estimation  for  finite  element  approximation  of  Stokes'  problem 

A  new  approach  to  a  posteriori  error  estimation  for  linear  elasticity  and  Stokes 
problems,  and  which  is  applicable  to  steady  state  Navier  Stokes  equations,  is  presented. 
The  approach  extends  earlier  work  |1)  to  problems  involving  an  incompressibility 
constraint.  The  method  is  valid  for  non-uniform,  irregular  and  adaptively  designed  h-p 
meshes  including  cases  where  the  polynomial  degree  need  not  be  uniform. 

The  approach  makes  use  of  duality  arguments  and  is  based  on  the  element  residual 
method  (ERM),  in  which  one  solves  small  local  problems  to  obtain  error  estimators  on 
each  element  in  the  mesh.  Key  featu,cs  of  the  approach  are  the  development  of  a 
systematic  approach  to  the  derivation  of  boundary  conditions  for  the  local  error  residual 
problems.  Significantly,  the  local  problem  need  not  in  fact  be  a  local  Stokes  type 
problem,  a  point  which  is  of  great  importance  in  simplifying  the  approximation  of  the 
local  problems. 

In  the  talk,  the  theoretical  foundations  of  the  method  are  outlined  and  results  of  its 
applications  to  several  representative  problems  are  presented. 


1. 


M  Ainsworth  and  J  T  Oden,  A  unified  approach  to  a  posteriori  error  estimation 
based  on  element  residual  methods.  Numeriscbe  Mathematik,  To  appear  (1992). 
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S  AM1NI  and  N  D  MA1NES 

Boundary  integral  equations  for  the  exterior  acoustic  problem 

The  boundary  integral  solution  of  the  Helmholtz  equation 

(V*  ♦  k1)*  -  0  .  p  €  D.  0) 

in  an  unbounded  domain,  D,,  exterior  to  a  smooth  C*  boundary  T,  satisfying  a  general 
boundary  condition  on  f  and  an  appropriate  radiation  condition  at  infinity  is  considered. 
By  a  careful  application  of  the  Green's  second  theorem  the  unique  solution  of  this 
problem  can  be  obtained  in  the  form 

*(P)  -  lrm~(PAXtrt  -  lrGk(p.q)^(qyrt,  p  e  D,.  (2) 

In  2-D,  the  Green’s  function  GJp,q)  =  i/4  Hjf'/k  |  p  -  q  |  ).  The  Helmholtz  integral 
representation  formula  (2)  requires  the  Cauchy  data  (<p,  d<p/Sn).  The  boundary  condition 
on  F  gives  either  <f>  or  d4>liin  or  a  relationship  between  the  two  in  the  case  of  "spring 
like"  scatterer.  We  need  another  boundary  relationship  between  <f>  and  d<t>/dn.  This  is 
usually  obtained  by  taking  the  limit  asp*  €  Dt  —  p  €  t  in  (2).  We  obtain  the  so-called 
"Surface  Helmholtz  Equation" 

-  (3) 

I  ^  ^ 

The  operators  on  either  side  of  equation  (3)  can  be  shown  to  be  singular  for  a  countable 
set  of  real  positive  values  of  k  which  correspond  to  the  standing  wave  solutions  of 
respective  interior  problems.  The  following  formulation  due  to  Burton  and  Miller  (1971) 
overcomes  this  non-uniqueness  problem  and  is  based  on  coupling  (3)  with  its  derivative 
in  the  direction  of  the  normal  to  T  at  p, 

( -~1  *  *  inN^fp)  -  (i,  ♦  pe  r.  (4) 

In  (4),  Lt  and  Mt  are  the  single  and  double  layer  potentials  with  M‘k  and  Nk  as  their 
respective  derivatives.  In  appropriate  Sobolev  spaces  hr(T)  it  can  be  shown  that  M„ 
M'k :  hr  -*  hr'1  are  smoothing  operators  of  order  -l  whilst  Nk :  hr "  -*  hT  is  principally 
an  order  -t- 1  differentiation  operator. 

In  this  paper  we  study  the  spectral  properties  of  the  four  Helmholtz  integral  operators. 
Ihe  choice  of  the  coupling  parameter  rj  and  also  the  convergence  analysis  of  collocation 
and  discrete  collocation  solutions  of  the  strongly  elliptic  pseudo-differential  operator 
equation  (4)  will  be  discussed.  We  also  report  on  several  2-grid  iterative  methods  for 
our  boundary  element  equations. 
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B  ANDERSSON 

On  the  modelling  errors  in  fracture  analysis  of  delaminated  plates 

Composite  materials  are  frequently  used  in  engineering  applications.  One  critical  failure 
mechanism  for  this  class  of  materials  is  related  to  the  delamination  of  sets  of  laminate 
layers.  For  safe  design  using  these  materials,  reliable  methods  for  prediction  of  growth 
and  arrest  of  delamination  damage  are  needed.  Delaminations  have  been  observed  to 
propagate  in  very  complex  patterns.  Initially  circular  delamination  fronts  may  under 
compressive  loading  (so-called  buckling-driven  delamination)  form  "worm’-like  shapes 
or  "telephone  cord’-shapes.  The  numerical  models  that  attempt  to  simulate  this  type  of 
delamination  growth  are  generally  based  on  nonlinear  plate  theory  and  fracture 
mechanics.  The  reasons  for  using  plate  models  are  that  the  dctaminations  appear  to  be 
"thin"  and  that  a  three-dimensional  analysis  is  considered  <o  be  too  costly.  A  critical 
question  in  this  context  is,  under  what  circumstances  do  plate  models  yield  accurate 
solutions? 

The  growth  of  delamination  fronts  is  assumed  to  be  governed  by  a  local  growth 
criterion  of  the  type  f(K/y),  K,/y),  K„Jy))  =  0,  where  Km(y)  are  the  stress  intensity 
functions  andy  is  a  coordinate  along  the  delamination  front.  It  is  of  special  importance 
to  separate  the  individual  fracture  modes  a,  since  experiments  indicate  significant 
differences  in  energy  release  rates  in  opening  and  sliding/tearing  modes.  The  reliability 
of  solutions  obtained  from  this  type  of  numerical  model  depend  critically  on  the 
modelling  error,  the  discretisation  error,  and  the  extraction  procedures  used  to  derive 
the  stress  intensity  functions  KJy).  For  plate  models,  bending  moments  and  shear  forces 
are  determined.  With  known  moments  and  forces,  the  stress  intensity  functions  are 
estimated  using  analytical  solutions  (split-beam  solutions)  for  a  two-dimensional  domain 
representing  the  delamination-tip-region. 

In  the  present  paper,  the  modelling  errors  obtained  using  the  Kirchoff  and  the 
Reissner-Mindlin  plate  models  are  studied.  Laminates  of  isotropic,  bi-material  isotropic 
and  orthotropic  material  having  quadratic,  elliptic,  or  (some)  real-life  shaped 
delamination  fronts  are  considered.  Three-dimensional  almost  exact  reference  solutions 
(i.e.  maximum  relative  pointwise  errors  in  K.(y)  of  the  order  10  ’)  are  derived  using  the 
h-p  version  of  FEM  and  a  novel  computational  procedure  for  extracting  edge  stress 
intensity  functions  at  3D  interfacial  delaminations  in  ortbotropic  laminates.  The 
solutions  for  the  Kirchoff  and  the  Reissner-Mindlin  plate  models  are  derived  using  the 
p  version  of  FEM.  Modelling  errors,  thus  isolated  from  discretization  errors,  are 
determined  virtually  exactly  for  different  plate  models,  fracture  mechanics  models, 
material  combinations,  delamination  shapes,  and  delamination/thicknesses  ratios  for  a 
simple  loading  case. 

Since  boundary-  layers  are  very  different  for  different  models,  so  are  calculated 
intensity  functions.  For  delaminations  with  a  diameter/thickness  ratio  of  *  =  20  and 
simple  shapes  of  the  delamination  front  (i.e.  quadratic  or  elliptic  shape),  the  modelling 
errors  (K’,D  -  K^“)IK’,°  are  of  the  order  5-15%.  For  near-circular  delaminations  with  k 
=  20  and  a  ’wavy"  boundary  (as  sometimes  observed  in  practice),  the  modelling  error 
may  be  as  large  as  25-40%.  The  conclusion  is  that  since  energy  release  rates  G,  are 
proportional  to  plate  model  solutions  often  do  not  provide  reliable  quantitative 
results  when  applied  to  real-life  delaminations  with  complex  boundaries. 
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P  ANDERSON,  J  C  ELLIOTT  and  S  E  P  DOWKER 

A  finite  difference  method  of  determination  of  local  coefficients 
in  dental  enamel  and  other  inhomogeneous  permeable  media 

Methods  to  study  transport  processes  in  teeth  (composed  of  inhomogeneous,  permeable 
materials)  are  important  in  the  understanding  of  carious,  destructive  and  repair 
processes.  In  normal  and  diseased  enamel,  and  many  other  inhomogeneous  materials, 
the  local  porosity  and  tortuosity,  and  therefore  the  diffusion  coefficient  vary  with 
position.  To  measure  diffusion  coefficients  in  homogeneous  permeable  materials  (e.g. 
a  glass  frit),  a  solute  with  high  X-ray  absorbance  (e.g.  1  mole  L 1  KJ)  is  allowed  to  diffuse 
from  a  constant  source  at  one  surface  into  water  contained  within  its  pores.  The 
concentration  of  Kl  at  any  point  along  the  1-dimensional  diffusion  path  can  be  calculated 
from  X-ray  attenuation  measurements.  Repeated  measurements  along  the  diffusion  path 
yield  a  time  series  of  concentration-distance  profiles  from  which  the  diffusion  coefficient 
for  KI  can  be  calculated  using  the  Crank-Nicolson  (CN)  algorithm  for  Fick’s  2nd 
equation,  the  aim  of  this  study  is  to  extend  this  method  to  measure  the  local  diffusion 
coefficient  VD,„,)  in  inhomogeneous  materials.  Experimental  data  was  simulated  for  the 
diffusion  of  a  solute  within  a  permeable  material  of  predetermined  inhomogeneity  using 
the  CN  algorithm  with  very  small  distance  and  time  steps.  Using  this  data,  D^,  was 
calculated  at  many  positions  using  a  modification  of  the  CN  algorithm.  Systematic  errors 
due  to  the  approximations  required  for  the  finite  difference  method  were  investigated. 
Errors  associated  with  intrinsic  uncertainties  due  to  X-ray  photon  counting  statistics  were 
studied  using  Monte-Carlo  simulations.  Improvement  in  accuracy  by  increasing  count¬ 
time  is  constrained  by  the  change  in  concentration  with  time.  However,  KI  can  be 
repeatedly  diffused  in  and  out  of  the  permeable  solid  and  an  averaging  algorithm  used 
which  effectively  increases  the  count-time  for  any  particular  concentration  determination 
without  introducing  systematic  errors. 


T  APEL 

Anisotropic  finite  elements  and  application  to  edge  singularities 

The  classical  local  interpolation  enor  estimates  (see  e.g.  Ciarlet  1978)  were  derived 
under  an  assumption  which  is  in  two  dimensions  known  as  Zlimal’s  minimal  angle 
condition.  This  condition  was  weakened  by  different  authors  (Jamet  1976,  Babuika/Aziz 
1976,  Kfi£ek  1989)  to  a  maximal  angle  condition.  But  the  possible  advantage  of  using 
mesh  sizes  with  different  asymptotics  in  different  directions  which  leads  to  small  angles, 
was  not  extracted. 

In  this  paper,  anisotropic  interpolation  error  estimates  in  two  and  three  dimensions, 
proved  by  Dobrowolski  and  Apel  (1992)  are  reviewed.  Here,  one  derives  benefit  from 
the  different  asymptotics  of  the  mesh  sizes. 

The  results  are  applied  to  the  finite  element  approximation  of  elliptic  equations  on 
domains  with  edges.  Practical  computations  show  that  the  proposed  finite  element 
meshes  are  well  suited  for  elliptic  problems  with  solution  of  r‘-type. 
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D  N  ARNOLD 

The  boundary  layer  for  the  Reissner-Mindlin  plate  model 

Tbe  solution  to  tbe  Reissner-Mindlin  plate  model,  in  contrast  to  that  of  the  bihannonic 
model,  exhibits  a  complex  dependence  on  the  thickness  of  the  plate.  For  thin  plates 
there  may  be  a  boundary  layer,  the  existence  and  structure  of  which  depends  on  the 
boundary  conditions,  the  plate  geometry,  and  the  solution  compbnent  considered.  For 
example,  the  transverse  displacement  variable  does  not  exhibit  any  edge  effect,  but  the 
rotation  vector  exhibits  a  boundary  layer  except  for  very  special  data  or  plate  geometry. 
The  bending  moment  tensor  and  shear  force  vector  have  more  pronounced  boundary 
layers.  The  boundary  layer  is  stronger  for  a  free  or  simply-supported  plate  than  for  a 
clamped  one.  In  this  talk  we  describe  recent  results  which  furnish  a  complete  asymptotic 
description  of  the  dependence  of  the  solution  on  plate  thickness.  The  techniques  used 
to  obtain  the  asymptotic  expansions  and  rigorous  error  bounds  for  them  are  described 
and  illustrated  by  concrete  examples.  Applications  to  finite  element  analysis  are 
presented. 


A  AUGE  and  G  LUBE 

A  stablized  Galerkin  finite  element  method  for  solving  the  Boussinesq  equations 

We  consider  the  Boussinesq  approximation  of  tbe  stationary,  incompressible  and  n on- 
isothermal  Navier-Stokes  equations.  Spurious  numerical  solutions  of  standard  mixed 
Galerkin  finite  element  formulations  may  be  generated  by 

(a)  inappropriate  combinations  of  velocity/pressure  interpolation  functions  and/or 

(b)  the  presence  of  (locally)  dominating  convective  terms. 

As  a  remedy  we  extend  the  least-squares  stabilization  approach  [FF]  to  the  Boussinesq 
equations.  In  particular,  we  add  weighted  residuals  of  the  basic  equations  to  the 
Galerkin  discretization  in  order  to  stabilize  tbe  method  and  to  allow  arbitrary  finite 
element  approximations. 

For  a  linearized  problem  of  Oseen  type,  we  prove  asymptotic  error  estimates  on  a 
non-uniform  mesh  and  give  an  explanation  of  the  parameter  design  which  is  frequently 
used  in  CFD.  As  a  basic  tool  we  take  advantage  of  a  modified  Babuska-Brezzi 
condition.  Furthermore,  we  extend  existence  and  convergence  results  to  the  case  of 
regular  solutions  of  the  original  problem  using  Banach’s  fixed  point  theorem  |LA) 

We  conclude  with  numerical  results  for  some  benchmark  problems. 

(FF)  Franca,  L  P,  Frey,  S  L,  Hughes,  T  J  R:  Comp.  Meths.  Appl.  Mech.  Engrg.  95  S 
278  (1992). 

(LA)  Lube,  G,  Auge,  A:  Regularized  mixed  finite  element  approximations  of 
incompressible  flow  problems.  Preprint  TU  Magdeburg  1992. 
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S  AZIMI 

A  new  elemental  formulation  for  the  dynamic  analysis  of  structures 

Use  of  a  new  "Dynamic  Mode  Shape  Function"  (DMSF)  based  on  the  dynamic  behaviour 
of  every  structural  element  is  under  study.  Using  this  DMSF  the  new  elemental  stiffness 
and  mass  matrices  are  obtained.  The  solution  of  the  problem  for  the  free  vibration 
results  in  exact  frequencies  and  mode  shapes  of  the  structures.  As  an  application  the 
method  has  been  applied  to  a  bar  finite  element  Some  examples  of  vibration  of  bar 
systems  have  been  considered  for  clarification  and  the  results  are  compared  with  the 
previous  published  results. 


C  BAILEY,  M  CROSS  and  P  CHOW 

A  finite  volume  procedure  to  predict  deformation  and  residual  stress  in  castings 

The  foundry/casting  industry  plays  a  major  part  within  the  manufacturing  base  of  the 
United  Kingdom.  For  a  new  cast  design  the  present  practice  is  to  undertake  numerous 
trials  and  tests  to  find  the  optimum  design  options.  This  practice  is  very  wasteful  in  both 
materials  and  lost  man  hours. 

During  the  casting  process  hot  liquid  metal  is  poured  into  a  mould.  Once  the  mould 
is  filled  with  liquid  metal,  the  casting  together  with  the  running  system  and  feeders  will 
cool  and  solidify.  The  development  of  a  code  to  completely  simulate  all  the  macroscopic 
processes  involved  in  castings  will  greatly  aid  the  fou  ndry  engineer  in  reducing  lead  times 
and  producing  sound  castings.  The  developed  code  must  include  the  following  types  of 
analysis: 

1)  Fluid  flow  analysis  for  mould  filling. 

2)  Thermal  flow  analysis  of  solidification  and  the  evolution  of  latent  heat 

3)  Fluid  Flow  analysis  to  represent  thermal  convection. 

4)  Deformation,  hence  stress  analysis  to  predict  geometrical  soundnessof  cast. 

5)  Macroscopic  porosity  formation. 

There  is  a  high  degree  of  coupling  between  the  above  phenomena,  for  example  the 
convective  currents  are  dependent  on  the  temperature  changes.  These  in  turn  are 
dependent  on  the  geometric  deformation  as  beat  losses  are  restricted  across  the 
cast/mould  interface  due  to  gap  formation.  Obviously  it  would  be  of  great  benefit  to 
solve  all  the  above  Jsing  a  single  code. 

The  route  taken  in  this  research  is  to  discretise  the  conservation  equations, 
representing  the  above  phenomena,  using  Finite  Volume  procedures.  This  presentaiton 
will  concentrate  on  the  methods  used  for  deformation  and  stress  and  how  these  methods 
have  been  coupled  to  the  flow  and  solidification  calculations  to  simulate  some  simple 
castings. 


M  I  BAINES 


Adaptive  methods  for  grid  generation  and  finite  element  solutions 
of  time-dependent  partial  differential  equations 

A  direct  variational  approach  has  been  used  to  generate  algorithms  which  determine 
optimal  discontinuous  piecewise  linear  and  constant  L,  fits  with  adjustable  nodes  to  a 
continuous  function  (1).  The  process  generates  data-dependent  grids.  The  algorithms 
are  fast  and  robust,  the  mesh  cannot  tangle,  and  convergence  can  be  proved  under 
certain  conditions. 

It  can  be  shown  that  one  iterate  of  such  an  algorithm  corresponds  to  an  explicit  time 
step  of  the  Moving  Finite  Element  (MFE)  procedure  for  the  solution  of  a  particular 
time-dependent  differential  equation  [2].  When  the  MFE  procedure  is  run  to  steady 
state  for  this  equation  it  follows  a  path  leading  to  the  best  fit,  with  adjustable  nodes,  erf 
the  steady  state  operator  to  the  zero  function.  The  result  generalises  to  a  wider  class  of 
differential  equations  and  operators. 

For  first  order  partial  differential  equations  it  is  known  [3]  that  the  MFE  procedure 
approximates  the  method  of  characteristic  strips.  The  above  result  shows  that  there  also 
exists  a  superimposed  speed  component  arising  from  the  L,  projection  which  is  "best  fit 
seeking". 

An  alternative  moving  grid  strategy  is  proposed  based  on  a  two  pass  method  [4].  The 
first  pass  is  a  fixed  grid  predictor:  this  is  followed  by  one  or  two  iterations  of  the  best 
fit  algorithm  with  adjustable  nodes.  A  second  pass  of  the  solver  is  then  made  on  the 
moving  grid,  the  motion  being  provided  by  the  above  node  adjustment 
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A  BALOCH,  P  TOWNSEND  and  M  F  WEBSTER 
Viscoelastic  effects  in  expansion  flows 

In  this  work  we  investigate  expansion  flows  in  both  two-dimensional  and 
three-dimensional  domains.  We  are  particularly  interested  in  comparing  numerical 
solutions  for  Newtonian  and  viscoelastic  fluids  over  a  range  of  inertial  and  elasticity 
values.  A  comparison  is  made  of  our  numerical  solutions  with  experimental  results,  and 
some  conclusions  are  reached  with  respect  to  vortex  enhancement,  and  presence  and 
origin  of  vortices.  For  the  two-dimensional  study  a  40:3  ratio  geometry  is  considered  in 
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line  with  the  equivalent  experiments  and  solutions  are  tested  for  accuracy  and 
convergence  on  more  than  one  mesh.  In  the  three-dimensional  case  the  40:3  expansion 
ratio  is  extended  into  the  third  dimension  to  provide  a  40:3:3  geometry  and  both 
Newtonian  and  viscoelastic  solutions  are  computed. 

The  numerical  scheme  is  a  time-stepping  Taylor-Galerltin  finite  element  method  that 
involves  fractional  stages  at  each  time  step  through  the  additional  consideration  of  a 
pressure-correction  scheme.  Such  an  algorithmic  approach  is  implemented  in  a 
semi-implicit  fashion,  and  is  accurate  and  robust  when  solving  Newtonian  flows. 
However,  for  the  viscoelastic  regime  more  sophistication  is  required  to  overcome 
nonlinear  instabilities  that  arise.  This  necessitates  the  further  consideration  of  a 
consistent  streamline  upwind  weighted  residual  approach.  Also  the  constitutive  models 
are  selected  from  a  generalised  class  of  Pban-Thien/Tanner  fluids  which  exhibit 
shear-tbinningand  a  finite  extensions!  viscosity.  Comments  are  made  as  to  the  suitability 
of  these  particular  choices. 


A  BALOCH,  P  TOWNSEND  and  M  F  WEBSTER 

Influence  of  Extensional  Viscosity  on  Vortex  Development  in  Complex  Geometries 

In  this  study  we  are  concerned  with  the  influence  of  strain  thickening  and  inertia  upon 
the  steady  laminar  incompressible  flow  of  non-Newtonian  fluids.  We  concentrate  on 
flows  through  circular  contraction  geometries  and  chose  to  work  at  contraction  ratios  of 
four  to  one.  A  generalised  Newtonian  material  description  is  adopted  with  an  extensional 
viscosity  depending  on  second  and  third  invariants  of  the  rate  of  deformation  tensor. 
Particular'attention  is  paid  to  the  effect  of  the  abrupt  reentrant  corner  by  introducing 
in  contrast  a  rounded  corner  geometry.  Under  such  circumstances  this  work  casts  light 
on  the  formation  of  vortices,  their  origin  and  subsequent  growth  with  increase  in  flow 
rate  and  elasticity  number.  Both  salient  and  lip  vortices  are  observed  under  different 
conditions. 

The  numerical  scheme  employed  in  the  simulations  is  a  semi-implicit  time-stepping 
procedure  that  combines  the  merits  of  a  Taylor-Galerkin  weighted  residual  method  with 
those  of  a  pressure-correction  method.  The  resulting  procedure  is  a  fractional  staged 
scheme  taken  over  each  time  step,  that  has  already  proved  its  flexibility  in  solving  a 
range  of  complex  incompressible  viscous  flows.  The  required  extensional  behaviour  is 
incorporated  through  an  appropriate  functional  form  for  the  viscosity.  Our  numerical 
solutions  are  compared  with  other  numerical  predictions  in  the  literature,  and  with  some 
experimental  data. 


M  BERNADOU 

On  the  approximation  of  junctions  between  thin  shells 
by  curved  C' -finite  element  methods 

Many  thin  plate  problems  are  formulated  on  curved  boundaty  plane  domains.  Likewise, 
by  using  a  mapping  technique,  thin  shell  problems  can  be  formulated  on  such  curved 
boundaty  plane  domains.  Such  mapping  techniques  can  be  extended  to  the  junction 
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between  thin  shells.  This  means  that  for  all  these  cases,  we  need  to  approximate  a 
fourth  order  problem  set  on  curved  boundary  plane  domains.  These  approximations  by 
conforming  finite  element  methods  require  the  use  of  curved  C, -finite  elements. 

In  this  talk,  we  will  present  some  results  obtained  in  the  analysis  of  curved  C'-finite 
elements  and  on  their  applications  to  solve  the  above  mentioned  problems. 
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M  BERZINS  and  J  WARE 

Reliable  finite  volume  methods  for  time-dependent  PDEs 

In  the  numerical  solution  of  time-dependent  pde's  the  accuracy  is  influenced  by  the 
spatial  discretisation  method  used,  the  spatial  mesh  and  the  method  of  time  integration. 
In  particular  the  spatial  discretisation  method  and  position  of  the  spatial  mesh  points 
should  ensure  that  the  spatial  error  is  controlled  to  meet  the  user’s  requirements.  It  is 
then  desirable  to  integrate  the  o.d.e.  system  in  time  with  sufficient  accuracy  so  that  the 
temporal  error  does  not  corrupt  the  spatial  accuracy  or  the  reliability  of  the  spatial  error 
estimates.  However,  in  most  existing  methods  for  time  dependent  p.d.e.’s  either  a 
stability  control  and  a  specialised  time  integration  method  is  employed  or  an  o.d.e.  solver 
is  used  which  makes  use  of  local  time  error  control  that  is  unrelated  to  the  spatial  error. 

This  paper  is  concerned  with  combining  the  spatial  and  temporal  error  balancing 
approach  of  Berzins,  [5],  with  the  adaptive  mesh  algorithms  of  Berzins  et  al.  [1],  [3]. 
TTiis  combination  of  error  control  strategies  provides  an  algorithm  that  adapts  the  space 
mesh  in  order  that  the  spatial  error  is  controlled  and  then,  accordingly,  adjusts  the  local 
error  tolerance  used  in  the  o.d.e.  integrator  so  that  the  two  errors  are  related  in  some 
way.  This  not  only  ensures  that  the  main  error  present  is  that  due  to  spatial 
discretisation  but  offers  the  tantalising  possibility  of  overall  error  control. 

In  the  error  balancing  procedure  the  local  time  error  is  controlled  so  that  it  is  a 
fraction  of  the  growth  in  the  spatial  error  over  the  timestep.  An  analysis  is  presented 
to  show  that  although  the  method  attempts  to  control  accuracy  a  Courant-like  stability 
condition  is  also  satisfied.  Numerical  experiments  show  that  the  error  control  strategy 
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appears  to  offer  an  effective  method  of  balancing  the  spatial  and  temporal  errors. 

The  unstructured  triangular  mesh  spatial  discretisation  method,  see  [1]  and  [4],  is  a 
cell-centered,  finite  volume  scheme  that  is  shown  to  achieve  second-order  accuracy  by 
use  of  a  six  triangle  stencil  around  each  edge.  The  growth  of  the  spatial  error  is 
measured  locally  in  time  by  using  the  computed  solution  as  the  high-order  solution  and 
constructing  an  error  estimate  by  calculating  a  computationally  inexpensive  low-order 
solution.  Thus  local  extrapolation  in  space  is  effectively  being  used.  These  spatial  error 
estimates  are  shown  experimentally  to  mimic  the  behaviour  of  the  global  spatial  error. 
Adaptive  unstructured  mesh  techniques  e.g.  [1]  or[3],  are  used  to  rezone  the  mesh  at 
times  chosen  by  the  time  integration  algorithm.  Selection  of  appropriate  times  is  made 
by  using  a  combination  of  present  estimated  errors  and  predicted  future  errors.  The 
time  integration  technique  used  in  this  paper  is  a  new  adaptive  Theta  Method,  [2], 
suitably  modified  to  use  the  new  error  control  approach. 

The  implementation  of  the  algorithm  raises  a  number  of  issues.  For  instance,  since 
the  accuracy  tolerance  for  the  time  integration  over  the  next  time  step  after  spatial 
remeshing  depends  partially  on  the  error  incurred  prior  to  spatial  remeshing,  this 
tolerance  must  be  modified  according  to  the  expected  reduction  in  the  spatial 
discretisation  error. 

The  present  algorithm  thus  goes  some  way  towards  providing  a  reliable  solution 
method  for  two-dimensional  time-dependent  flows.  The  numerical  results  show  that  this 
can  be  achieved  for  both  simple  convection  type  problems  and  for  Euler  flows  in  two 
space  dimensions. 
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I  BOCK  and  J  LOVISEK 

On  the  discretization  of  an  optimal  design  problem  for 
an  unilaterally  supported  viscoelastic  plate 

The  discretizaiton  of  an  optimal  design  problem  for  a  visco-elastic  plate  with  respect  to 
its  variable  thickness  is  solved.  The  middle  surface  of  the  plate  is  identified  with  an 
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open  bounded  plane  domain  C  with  a  Lipschitz  boundary  dG  -  T.  U  iy  The  plate  is 
clamped  on  T,  and  unilaterally  supported  on  iy  Let  K  C  HP(Cl)  be  the  closed  convex 
set  of  admissible  deflections.  The  deflection  w(e,Ux)  is  a  solution  of  the  initial-boundary 


value  problem  for  a  pseudoparabolic  variational  inequality. 
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We  connect  with  (1),  (2),  (3)  the  design  problem 
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where  U*  C  Cf'(G)  is  the  admissible  set  of  thickness-functions,  X  is  any  Hilbert  space, 
D  :  X  -*  H2(G)  is  a  linear  bounded  operator.  Let  K*  and  Iff  be  finite  element 
approximations  of  the  sets  K  and  Ua  respectively,  d  =  TIN.  Using  finite  differences  with 
the  step  d  with  respect  to  t  we  approximate  the  problem  (1)  -  (4)  by  N  elliptic  optimal 
control  problems  with  admissible  deflections  and  thickness-functions  from  Kk  and  Uu 
respectively. 


L  BRUSA  and  F  RICCIO 

A  hybrid  direct-iterative  method  for  the  solution  of  finite  element 
linear  systems  with  positive  definite  matrices 

The  major  cost  of  many  finite  element  analyses  arises  from  the  solution  of  large  sparse 
sets  of  linear  equations  which  can  often  be  performed  only  by  using  vector  and  parallel 
computers.  Direct  methods,  based  on  Gaussian  elimination,  are  efficient  but  they  can 
easily  require  prohibitively  large  amounts  of  storage,  which  hamper  their  practical 
application  for  complex  engineering  analyses.  This  drawback  is  avoided  by  using 
iterative  schemes,  mainly  the  conjugate  gradient  method,  but  their  performance  is  often 
unsatisfactoiy  because  of  the  large  number  of  iterations  required  for  convergence,  due 
to  the  ill-conditioning  of  the  matrices. 

Hybrid  methods,  based  on  the  use  both  of  direct  and  iterative  techniques,  have  been 
proposed  in  an  attempt  to  exploit  the  advantages  of  the  two  approaches  and  to  mitigate 
their  drawbacks.  The  starting  point  of  the  methods  considered  in  the  paper  is  the 
substructure  technique  which  arranges  the  elements  in  groups  (or  substructures)  so  that 
two  classes  of  variables  can  be  defined,  i.e.  the  internal  veriables  relevant  to  nodes  within 
substructures,  and  the  interface  variables  relevant  to  nodes  belonging  to  two  or  more 
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substructures. 

The  internal  variables  are  decoupled  and  are  eliminated  by  applying  a  direct  method 
to  factorize  the  related  matrices.  The  reduced  system  obtained  after  this  step  involves 
only  the  interface  variables  and  is  solved  by  means  of  the  preconditioned  conjugate 
gradient  method  (PCGM).  In  this  way  the  efficiency  of  Gaussian  elimination  is  exploited 
in  the  solution  of  medium-sized  problems,  while  PCGM,  applied  to  solve  the  reduced 
system,  deals  with  a  matrix  which  is  better  conditioned  than  the  original  one.  This 
feature  of  the  reduced  system,  as  far  as  we  know,  has  been  experimentally  verified  and 
demonstrated  only  for  special  matrices. 

The  paper  presents  a  proof  that  the  condition  number  of  the  Schur  complement  (the 
matrix  of  the  reduced  system)  of  a  positive  definite  matrix  is  lower  than  the  one  relevant 
to  the  global  matrix.  On  the  basis  of  this  result,  indications  are  given  to  build 
preconditioners  for  the  reduced  system,  starting  from  those  defined  for  the  original 
problem.  Some  of  these  preconditioners  are  tested  for  the  solution  of  structural 
mechanics  problems  and  the  performances  are  compared  with  those  obtained  by  applying 
PCGM  to  the  complete  problem.  The  numerical  experimentation  has  been  performed 
on  the  parallel  shared  memory  multiprocessor  ALLIAN  FX/80  and  future  work  will 
concern  the  implementation  of  the  method  on  a  distributed  memory  hardware  platform. 


T  BULENDA 

Hybrid  Amoldi- Newton  Algorithm 

Applying  a  predictor-corrector  scheme  to  the  augmented  set  of  nonlinear  equilibrium 
equations  in  statics  G(u,A)  =  0,  we  get  for  an  inexact  Newton-Raphson  iteration  in  the 
corrector  phase  (superscripts  stand  for  the  iteration  number): 

* 

(•> 

with  the  tangential  stiffness  matrix  'KT  =  '(dR/du),  the  residuum  vector  R,  the  vector  of 
inner  forces  F  (dependent  on  the  vector  of  displacements  u),  the  reference  load  vector 
P  and  the  scaling  factor  for  the  loads  X.  The  equation  f  =  0  describes  the  surface  of 
iteration  and  f„  =  \dfjdu),  %  =  (d(JdX).  The  inexact  Newton  method  requires  some 
coupling  of  the  linear  equation  solver  and  the  Newton  iteration.  This  is  achieved  via  the 
term  '8  which  sets  the  residuum  *r  of  the  linear  equation  in  relation  to  that  of  the 
Newton  process  G.  The  value  of  “8  is  chosen  as  a  forcing  sequence.  In  the  present  work 
we  use:  ‘4  =  min  (1  O’,  0.1  |j  G  jj  /U). 

To  get  an  approximate  solution  of  the  linear  equations,  we  need  a  proper  equation 
solver.  In  this  work  the  Amoldi  algorithm  is  used.  It  can  be  regarded  as  a 
generalization  of  the  well  known  Lanczos  algorithm  to  unsymmetric  system  matrices. 
The  solution  of  our  linear  system  is  thus  reduced  to  the  repeated  computation  of 
(system-matrix)-vector  products.  Via  preconditioning  and  reorthogonalization 
respectively  the  algorithm's  performance  and  stability  is  improved.  As  the  Amoldi 


•r  -  -  F(u)*‘kr 


-AX'A) 


+  ir  *  -  Q  +  lr  with  Jgj  &  it 


19 


algorithm  can  handle  unsymmetric  matrices,  we  solve  directly  the  augmented  system  (1) 
instead  of  using  the  widely  known  partitioning  method  of  Batoz  and  Dhatt  (1979). 

The  results  of  our  computations  show  great  savings  in  the  total  number  of  Arnoidi 
iterations  we  get  by  working  with  the  augmented  system  instead  of  the  partitioning 
method.  Whether  this  is  a  saving  in  computer  time  also,  depends  on  the  structure  of  Kj 
and  the  cost  relation  for  building  the  matrix-vector  product  in  each  step  of  the  Arnoidi 
process  and  performing  one  Arnoidi  step.  If  KT  is  unsymmetric,  the  augmented  system 
is  favourable  anyway.  Then  the  additional  effort  for  solving  an  n+ 1  -dimensional  system 
can  be  neglected  compared  with  solving  a  second  n-dimensional  system  when  using  the 
partitioning  method.  If  KT  is  symmetric,  the  decision,  which  way  to  choose,  depends  on 
the  above  mentioned  cost  relation.  The  higher  the  cost  for  the  matrix-vector  product, 
the  more  favourable  is  the  augmented  system.  For  any  strategy  the  inexact  Newton 
method  additionally  produces  great  savings. 


M  BUONSANTI 

A  FEM  Approach  for  the  Analysis  of  Composite  Elements 

The  study  of  the  behaviour  under  load  of  heterogeneously  built  solids  is  treated  by 
analysing  the  contact  area  of  two  different  constituent  details.  The  problem  of  unilateral 
contact  with  or  without  friction  as  well  as  the  deformative  adaptation  that  results  from 
the  application  of  the  loads  is  studied.  The  disjunction  among  elements,  both  in  the 
particular  and  bi-dimensional  case  is  probed.  In  each  case  we  suggest  a  FEM  model  for 
numerical  applications. 


P  BURDA 

Cubic  finite  elements  in  elliptic  problems  with  interfaces  and  singularities 

In  applications  of  the  finite  element  method,  essentially  two  approaches  are  most  widely 
used.  First  one  makes  use  of  some  kind  of  refinement  of  the  mesh  near  the  singularity. 
Second  one  uses  additional  trial  functions  in  the  neighbourhood  of  the  singularity.  In 
our  paper  we  give  an  alternative  to  these  techniques,  based  on  using  cubic  polynomials 
in  the  FEM. 


P  CARTRAUD,  C  WIELGOSZ  and  O  DEBORDES 
Computation  of  the  homogenized  behaviour  of  a  gasket 

The  present  study  is  devoted  to  the  determination  of  mechanical  behaviour  of  a  cylinder- 
head  gasket.  We  consider  here  a  gasket  made  of  an  heterogeneous  material,  with  an 
iron  core  (perforated  sheet)  and  a  paper  (non  asbestos  material)  which  is  clamped  by 
rolling  on  both  sides  of  the  sheet.  This  gasket  exhibits  an  elastoplastic  behaviour  and 
its  effective  properties  are  computed  by  an  bomogenzatkm  technique. 

The  homogenization  method  implies  a  good  description  of  the  microstnicture,  and 
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especially  the  characteristics  of  the  gasket's  components.  At  this  stage,  the  main 
difficulty  is  to  model  the  paper  behaviour.  This  behaviour,  displayed  with  compression 
tests,  looks  like  that  of  soil. 

We  therefore  use  an  elastoplastic  model,  including  non-linear  and  anisotropic 
elasticity,  plasticity  with  two  yield  surfaces  and  isotropic  hardening.  One  homogeneous 
test  is  used  to  determine  the  material  parameters  based  on  an  analytical  solution.  We 
then  compute  the  numerical  solution  of  another  compression  test,  with  the  finite  element 
code  SIC  (University  Technologique  de  Compi£gne),  which  incorporates  specific 
developments  for  the  numerical  integration  of  elastoplastic  constitutive  laws  for 
geomaterials.  We  obtain  good  agreement  between  the  model  and  experimental  results. 

One  can  then  compute  the  three  dimensional  homogenized  behaviour  of  the  gasket 
We  first  make  use  of  the  planar  periodicity  of  the  gasket  material  (due  to  the  sheet 
perforation)  to  reduce  the  study  to  that  of  a  representative  volume  element  (r.v.e.).  The 
homogenized  behaviour  of  the  gasket  is  then  obtained  from  numerical  simulation  of 
elementary  mechanical  tests  on  a  finite  element  model  of  this  r.v.e.  Here  again  we  use 
the  finite  element  code,  SIC,  in  which  a  specific  module  of  homogenization  has  been 
developed.  In  practice,  the  r.v.e.  is  submitted  to  an  average  strain,  increasing  in  time, 
and  we  compute  the  associated  average  stresses.  The  effective  properties  of  the  gasket 
define  the  relations  between  the  averages  of  strains  and  stresses.  One  can  thus  construct 
the  elastoplastic  evolution  law  of  the  gasket,  by  considering  successive  strains  in  different 
directions.  These  computations  compare  well  with  experimental  tests  on  the  gasket 

These  results  will  enable  us  to  incorporate  the  gasket  behaviour  more  effectively  into 
the  finite  element  computation  of  an  engine. 


I  C  CAVENDISH,  C  A  HALL  and  T  A  PORSCHING 

A  complementary  volume  approach  for  modelling  three-dimensional  Navier-Stokes  equations 
using  dual  Delaunay/Voronoi  tessellations 

In  this  paper  we  described  a  new  mathematical  approach  to  deriving  and  solving  co¬ 
volume  models  of  three-dimensional,  incompressible  Navier-Stokes  flow  equations.  The 
approach  integrates  three  technical  components  into  a  single  modelling  algorithm: 

1.  Automatic  Grid  Generation.  An  algorithm  previously  used  for  three-dimensional 
finite  element  mesh  generation  applications  is  used  to  automatically  discretize  the 
flow  domain  into  a  Delaunay  tetrahedronization  and  a  Voronoi  tessellation. 
These  two  geometric  constructs  provide  3-D  complementary  volume 
discretizations. 

2.  Co-volume  Equation  Generation.  Using  the  tetrahedronization  and  dual 
tessellation,  a  finite  volume  method  is  applied  to  derive  a  system  of  simple, 
second-order  convergent,  finite  difference  equations  to  model  velocity  and 
pressure  dependent  variables  satisfying  the  Navier-Stokes  equations. 

3.  Dual  Variable  Reduction.  A  graph  theoretic  analysts  technique  is  used  to 
transform  the  finite  difference  system  into  an  equivalent  banded  system  which  is 
considerably  smaller  than  the  finite  difference  system.  Unlike  most  Navier-Stokes 
solution  algorithms,  the  velocity  field  is  a  priori  discretely  divergence  free. 

We  display  a  nonsymmetric  flow  duct  of  rectangular  cross  section  in  which  flow  enters 
through  one  inlet  port  and  exits  through  two  outlets.  A  Delaunay  tetrahedronization  of 
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the  flow  domain  (boundary  faces  only)  into  1320  tetrahedra  is  produced  and  the  dual 
Voronoi  tessellation  which  contains  542  convex  Voronoi  polytopes  is  obtained.  The 
derived  covolume  difference  system  (which  we  call  the  primitive  system)  is  6003  x  6003 
while  the  dual  variable  system  is  approximately  one  third  as  large.  Numerical 
calculations  will  be  compared  with  flow  visualization  results  for  water  flow  through  the 
(transparent  plastic)  duct 
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P  CHAUSSECOURTE,  B  M&T1VET  and  G  NICOLAS 
Adaptive  mesh  refinement  for  the  3D  magnetic  code  TRIFOU 

For  some  years,  we  have  been  working  to  develop  TRIFOU,  a  finite  element  code 
devoted  to  3D  magnetic  field  problems.  One  application  is  the  modelling  of  eddy 
current  non  destructive  testing.  In  this  problem,  we  must  build  a  mesh  for  cracked 
bodies.  This  should  be  done  carefully:  an  accurate  representation  of  the  induced 
current  around  the  crack  is  necessary  to  achieve  a  good  representation  of  the  control 
signals.  But  forecasting  the  shape  of  die  currents  and  building  an  appropriate  mesh  are 
often  difficult.  Therefore,  we  decided  to  implement  an  adaptive  mesh  refinement 
strategy  in  TRIFOU. 

The  Maxwell  equations  for  time  periodic  situations  lead  to  solving  the  following 
probi'em:  find  the  complex  magnetic  field  b  satisfying  i  at  |i  b  +  curl  (curl(h)/o)  =  0, 
where  eo.fi  and  o  are  calculation  parameters.  The  variational  formulation  of  this 
problem  can  be  seen  as:  find  h  in  H  =  jh  £  1L*(E),  curl(h)  £  ILZ(E))  such  that 
a(h,h')  =  (  b,h'),  for  all  h'  in  H,  where  b  belongs  to  the  dual  of  H.  The 
approximate  solution,  hk,  is  sought  in  the  finite  dimensional  space,  Hk,  which  is 
connected  to  the  mesh,  Mk.  We  can  define  the  residual  of  the  discrete  formulation,  Rk, 
by:  (  Rk,h')  =  (  b,h')  -  a(hk,h'),  for  all  h'  in  H.  The  natural  norm  of  Rk  is  given 
by  the  maximum  of  |  <  Rk,b '  >  |  for  ail  b '  in  the  unit  sphere  of  H.  This  value 
is  a  norm  of  the  global  eiror,  but  its  computation  is  difficulL  Moreover,  this  value  does 
not  allow  us  to  say  where  and  how  the  mesh  should  be  refined.  Nevertheless,  we  shall 
use  this  idea  for  our  adaptive  mesh  refinement,  which  can  be  seen  as  an  optimisation 
algorithm:  the  successive  meshes  are  built  up  in  order  to  minimise  this  residual. 

The  mesh  Mk+I  is  obtained  by  refining  some  carefully  chosen  areas  of  the  mesh  Mk. 
We  define  a  unit  refinement  as  the  lowest  increase  in  the  approximation  space  Hk.  An 
adaptive  mesh  iteration  is  divided  into  two  stages.  First,  we  test  ail  the  possible  unit 
refinements  of  the  mesh:  for  each  of  these,  we  enrich  the  approximation  space  Hk 
temporarily  with  the  few  new  functions  related  to  this  virtual  refinement;  we  evaluate 
the  efficiency  of  this.  In  the  second  stage  of  the  process,  we  begin  by  selecting  the  most 
efficient  unit  refinements.  We  then  build  up  the  Hkt,  space,  adding  all  the  selected 
unit  refinements  to  Hk. 

To  evaluate  the  efficiency  of  a  virtual  unit  refinement,  we  consider  the  highest  value 
of  |  ( Rk,h ')  |  /  U  h '  |  for  h '  in  the  set  of  the  new  functions.  We  have  to  solve  a  small 
linear  system;  its  rank  is  approximately  the  number  of  degrees  of  freedom  which  are 
introduced  by  the  unit  refinement 

The  results  of  the  first  applications  of  this  technique  are  promising,  and  the 
refinement  occurs  in  the  immediate  surroundings  of  the  cracks  or  the  singularities. 
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11  CHEN 


Richardson  extrapolation  for  the  streamline  diffusion  finite  element  method 


The  Richardson  extrapolation  is  a  classical  technique  for  increasing  the  accuracy  in 
numerical  analysis  and  has  been  applied  to  the  standard  finite  element  method,  c.f. 
Blum,  Lin  and  Rannacher  (1].  The  purpose  of  this  paper  is  the  study  of  extrapolation 
methods  for  the  streamline  diffusion  finite  element  method  (the  SD  method)  that  has 
been  developed  by  Hughes  and  Brooks  [2]  for  solving  convection-dominated  convection- 
diffusion  problems  numerically.  The  theoretical  analysis  of  the  SD  method  has  been 
started  by  Johnson  and  Navert  in  (3).  Global  and  local  error  estimates  in  L2-norm  of 
order  0(hk*113)  have  been  derived.  These  estimates  are  optimal  for  general  quasi-uniform 
meshes.  In  Johnson,  Schatz  and  Wahlbin  [4]  the  pointwise  error  behaviour  was  analysed 
for  the  linear  finite  element  and  an  error  estimate  of  order  0(hiM  |  logA  |  w)  was  proved. 
An  improved  error  estimate  of  order  0(b“*  |  log/i  |  )  has  been  obtained  in  Niijima  (6). 
Up  to  now,  this  estimate  is  the  best  pointwise  error  estimate.  In  this  paper  we  derive 
error  asymptotic  expansions  with  remainders  of  order  up  to  0(hltm  |  logh  |  m)  for  the 
SD  method.  The^  are  the  theoretical  approach  of  the  extrapolation  methods  for 
increasing  accuracy.  C:ir  methods  require  that  the  exact  solution  is  sufficiently  locally 
smooth.  In  the  meantime,  a  quasi-optima!  pointwise  error  estimate  of  order 
0(h211 1  logh  |  ia)  will  be  derived  (or  locally  uniform  meshes  by  making  use  of  the  known 
sup»rconvergence  results  for  the  finite  element  interpolation. 

As  a  model  problem,  let  us  consider  the  following  boundary  value  problem. 


-  pn„  -  at*  *  uM  *  u  -  f  in  0 

*»l«  *  0  • 


(1) 


where  fl  e  R1  is  a  bounded  smooth  domain  and  p,  e  are  small  positive  parameters. 

Let  Sk  C  H'dfi)  be  a  piecewise  linear  finite  element  space  constructed  on  a  quasi¬ 
uniform  mesh  Tk.  The  so-called  streamline  diffusion  finite  element  approximation  u*  € 
Sk  to  the  problem  (I)  is  determined  through  the  relation 

-  (f.X  .  V  z  6  5,  .  (2) 


where 


B^u\x)  =  p(u,\x.)  ♦  c/n/.x,)  ♦  («*,*♦«*♦  P»X^ 

■  (Ps^PKa/.X,)  *  ♦  (1  -  P*X*,\x)  -» 

and  pk  =  ch,  ek  =  chm.  For  simplicity,  we  assume  that  e  s  ct,  p  s  pk. 
Then  we  will  derive  the  following  error  asymptotic  expansions 
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“4  “  *  +  «*«,  *  P*P^  *  k%  *  «*P*a,  *  «*&,  +  P*P£s 
+  ^Pu,  ♦  0(€+*J<W|lO*fc|W) 

Based  on  this  expansion,  we  design  some  extrapolation  methods  to  gain  higher  accuracy. 
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S  H  CHEN 

A  new  model  of  jointed  rock  masses  reinforced  by  passive,  fully-grouted  bolts 

Rock  bolts  have  been  used  as  general  reinfc  intent  measure  to  improve  the  strength 
of  jointed  rock  masses  in  recent  decades.  Several  useful  models  have  been  proposed  by 
some  researchers. 

This  paper  presents  a  new  model  for  the  jointed  rock  masses  reinforced  by  passive, 
fully-grouted  bolts.  The  model  is  based  on  an  "equivalent  continuum'  approach.  On  one 
hand,  the  increment  of  load  is  shared  among  the  bolts,  rock  materials,  and  joints 
respectively.  On  the  other  hand,  the  strain  increments  of  the  reinforced  jointed  rock 
masses  are  given  as  the  sum  of  the  increment  strains  of  the  reinforced  rock  materials  and 
reinforced  joints.  The  proposed  model  has  an  advantage  in  comparison  to  previous 
'equivalent  continuum"  models  -  it  can  calculate  the  bolt  stresses  separately  in  joints  and 
rock  materials.  Because  their  deformations  are  different,  this  is  a  more  reasonable 
approach  to  the  real  work  state  of  bolts  in  the  jointed  rock  masses.  A  computer 
program  of  the  Finite  Element  Method  has  been  compiled,  and  a  typical  case  study  of 
the  rock  slope  stability  problem  has  been  made. 
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Z  CHEN,  B  COCKBURN,  J  W  JEROME  and  C-W  SHU 


Runge-Kutta  local  projection  discontinuous  Galerkin  methods  for  conservation  laws 
and  applications  to  the  hydrodynamic  device  model 


In  this  paper  we  introduce  and  analyze  a  class  of  nonlinearly  stable  Runge-Kutta  local 
projection  discontinuous  Galerkin  (RKDG)  finite  element  methods  for  the  hyperbolic 
conservation  law 


u,  +  div  f(u)  =  0,  in  (0 ,T)  x  Cl 


(*I) 


with  suitable  initial-boundary  conditions,  where  Cl  C  ,  u  =  (u„  ....  um),  and  f  is  such 


d 

that  any  real  combination  of  the  Jacobian  matrices  'p  t  hasm  real  eigenvalues  and 

&  'an 


a  complete  set  of  eigenvectors.  The  main  difficulty  in  solving  (1.1)  is  that  solutions  may 
contain  discontinuities  even  if  the  initial-boundary  data  are  smooth.  Among  the 
successful  numerical  methods  for  solving  (1.1),  we  mention  the  nonoscillatoiy 
conservative  finite  difference  methods  such  as  total  variation  diminishing,  total  variation 
bounded,  and  essentially  nonoscillatory  schemes,  and  the  nonstandard  finite  element 
methods  such  as  streamline  diffusion  and  characteristic  Galerkin  approaches.  One 
distinctive  feature  of  the  methods  under  consideration  is  a  local  projection  limiting, 
which  comes  from  the  successful  nonoscillatory  finite  difference  methodology,  guarantees 
total  variation  boundedness  for  one-dimensional  nonlinear  scalar  equations  and  linear 
systems,  and  yields  a  local  maximum  principle  for  multi-dimensional  scalar  equations. 
The  main  advantages  of  these  methods,  over  most  finite  difference  methods,  are  that 
they  can  easily  handle  complicated  geometries  and  boundary  conditions,  and  that,  over 
most  other  finite  element  methods,  they  use  high  order  total  variation  diminishing 
Runge-Kutta  type  time  discretizations,  which  renders  them  explicit,  fully  parallelizable, 
and  computationally  efficient.  In  this  paper,  the  analysis  of  the  RKDG  methods  defined 
on  general  triangulations  is  carried  out  for  the  two-dimensional  system  (11),  d  =  2  and 
m  >  1.  It  is  shown  that  these  methods  satisfy  local  maximum  principles,  are  uniformly 
high-order  accurate,  and  converge  to  a  weak  solution. 

After  the  analysis  above,  the  application  of  the  RK.DGR  methods  to  the  hydrodynamic 
mode)  for  semiconductor  devices  is  presented: 


H,  ♦  div(nr)  •  0  , 

j,  <■  *  di v  p  ♦  p  .  V>  «  -«nK  -  *  (p,)t 

W,  *  div(rw)  «  «n*  .  E  -  div(r»7)  ♦  VT)  ♦  (w»)e  , 


where  n  is  the  electron  density,  r  is  the  velocity,  p  is  the  momentum  density,  t{>  0)  is 
the  electronic  charge,  E  is  the  electric  field,  T  is  the  temperaiuic.  w  is  (he  energy  density, 
k  is  the  heat  conduction  coefficient,  and  the  index  c  undoes  collision  terms.  It  is 
assumed  that  the  energy  bands  are  parabolic  as 
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3  1 

p  -  mnv,  w  m  ?.nT  +  .i/wrt|v|J  , 

where  m  is  the  effective  electron  mass.  The  electric  field  satisfies  the  Poisson’s  equation 
divfeV*)  -  -e(N0-NA-n),  B  -  -V* 

i 

where  N0  and  NA  are  the  densities  of  donors  and  acceptors,  respectively.  This  model 
defines  the  variables  n,  p,  w,  and  0,  treats  electron  flow  in  a  semiconductor  device  via 
the  Guler  equations  of  gas  dynamics,  with  the  addition  of  a  heat  conduction  term,  and 
is  a  coupled  system  of  hyperbolic,  parabolic,  and  elliptic  equations.  The  nonlinear 
hyperbolic  equations  support  shock  waves. 

Muted  finite  element  methods  are  used  for  the  approximation  of  the  electric  field. 
The  reason  for  including  the  mixed  methods  is  that  the  transport  equations  depend  on 
the  potential  only  through  this  field  and  these  methods  provide  better  approximation  of 
it  than  more  standard  Galerkin  approaches  would  give.  After  the  Lagrange  multipliers 
are  introduced  for  the  mixed  methods,  the  scheme,  a  combination  of  the  RKDG  and 
mixed  methods,  if  fully  parallelizable.  Extensive  numerical  simulations  have  been  carried 
out  on  the  Cray-2  at  the  Minnesota  Supercomputer  Center  and  the  Army  High 
Performance  Computing  Research  Center  in  USA. 


E  R  CHRISTENSEN 

Structural  dynamic  finite  element  analyses  of  space  shuttle  propulsion  components 

This  paper  describes  three  finite  element  structural  dynamic  analyses  of  the  space  shuttle 
main  engine  and  solid  rocket  boosters  which  have  recently  been  conducted  for  the 
NASA  Marshall  Space  Flight  Center  in  Huntsville,  Alabama,  USA,  in  support  of  the 
Advanced  Turbopump  Development  (ATD)  and  Advanced  Solid  Rocket  Motor  ( ASRM) 
programs.  These  analyses  were  performed  using  the  Engineering  Analysis  Language 
(EAL)  and  COSMIC  NASTRAN  finite  element  codes.  In  the  first  analysis,  three 
dimensional,  solid  finite  element  models  (FEM)  of  the  ATD  fuel  pump  and  LOX  pump 
first  stage  turbine  blades  were  constructed  and  analyzed  using  EAL.  Since  modal 
analysis  results  indicated  possible  resonance  conditions,  the  transient  response  of  the 
blades  due  to  the  operating  pressure  loads  was  calculated.  Several  transient  analyses 
were  conducted  for  various  operating  conditions  with  the  results  indicating  the  possibility 
of  fatigue  life  problems.  The  second  analysis  involved  an  investigation  of  the  fluid- 
structure  interaction  between  liquid  hydrogen  fuel  and  toroidal  shaped  inlet  structure  on 
the  high  pressure  fuel  turbopump.  A  FEM  representing  both  the  fluid  and  the  structure 
was  constructed  and  a  modal  analysis  was  performed.  The  results  indicated  the  presence 
of  a  large  number  of  "spurious”  fluid  modes  occurring  at  nonzero  frequencies.  These 
spurious  modes  are  a  direct  result  of  the  displacement  formulation  of  the  EAI  fluid 
elements  which  does  not  impose  the  constraint  of  irrotationality  on  the  fluid.  When 
these  constraints  were  enforced  through  the  use  of  a  global  penalty  function,  however, 
it  was  found  that  the  spurious  modes  could  be  removed  from  the  solution.  The  final 
analysis  involved  determining  the  effect  of  internal  pressure  on  the  dynamic 
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characteristics  of  the  ASTM  which  is  currently  being  developed  as  a  replacement  for  the 
solid  rocket  motors  presently  used  on  the  shuttle.  Several  new  NASTRAN  elements 
have  been  developed  to  model  the  pressure  stiffness  which  is  a  nonlinear  effect  due  to 
the  change  in  direction  of  the  pressure  force  as  the  structure  deforms.  Neglecting  this 
effect  can  lead  to  errors  of  up  to  48%  in  the  computed  frequency  of  the  pitch/roll  mode. 
Prediced  frequencies  obtained  using  the  new  NASTRAN  elements,  however,  agree  with 
modal  test  results  to  within  5%. 


R  K  COOMER  and  I  G  GRAHAM 

Domain  decomposition  techniques  for  massively  parallel  machines 

In  this  talk  we  describe  massively  parallel  domain  decomposition  methods  for  solving 
large  ill-conditioned  linear  systems  arising  from  discretisations  of  symmetric  elliptic 
PDEs.  The  method  has  particular  power  when  the  ill-conditioning  is  made  worse  by  the 
presence  of  discontinuous  coefficients  in  the  PDE.  Methods  based  on  mapping 
subdomains  to  an  array  of  processors  and  local  elimination  of  interior  unknowns  are 
reviewed.  The  resulting  Schur  complement  system  is  solved  by  the  preconditioned 
conjugate  gradient  method  with  preconditioner  (in  2D)  based  on  local  solves  in  vertex 
and/or  edge  spaces,  plus  a  coarse  grid  solve  to  simulate  global  interaction  of  nodes.  This 
is  an  "additive  Schwarz"  type  method,  and  its  analysis  has  been  recently  investigated  by 
M  Dryja,  O  Widlund  and  B  Smith.  It  is  explained  why  this  method  yields  an  optimal 
algorithm  (i.e.  its  convergence  rate  is  independent  of  the  number  of  degrees  of  freedom), 
but  in  general  its  convergence  rate  is  affected  by  the  jumps  of  the  coefficients  of  the 
PDE  acroks  subdomain  boundaries.  On  the  other  hand,  neglecting  the  vertex  solves 
yields  a  weakly  suboptimal  algorithm  but  with  convergence  rate  unaffected  by  these 
jumps.  The  method  is  implemented  for  a  class  of  model  problems  on  the  MasPar  MP-1, 
a  SIMD  machine  with  1024  processors  using  inner  iterations  to  solve  the  preconditioning 
problems.  Numerical  results  which  support  the  theoretical  properties  of  the  algorithms 
are  reported. 


A  CRAIG 

p-version  preconditioning  for  the  mass  matrix 

The  p-version  of  the  finite  element  method  has  become  more  popular  over  the  last  few 
years.  The  increase  in  popularity  is  due  to  several  factors;  including  the  greater 
return  in  accuracy  for  problems  with  singularities  on  the  mesh  and  the  case  of  design  of 
p- adaptive  codes.  Here,  instead  of  refining  the  mesh,  we  increase  the  accuracy  by  using 
higher  degree  basis  functions. 

There  is  an  observed  problem  which  causes  computational  difficulties:  the  discrete 
operator  from  the  finite  element  method  is  highly  ill-conditioned. 

Widlund,  Pasciak  et  aL  and  Babuika  et  aL,  to  name  but  a  few,  have  all  given 
experimental  results  which  reduce  the  order  of  the  condition  number  for  h  and 
p-versions  of  the  finite  element  stiffness  matrix  from  polynomial  to  log-squared  in  the 
number  of  degrees  of  freedom;  also,  numerical  results  suggest  a  similar  method  for  the 
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hybrid  h-p-vc rsion. 

An  additional  problem  arises  in  the  p- version.  We  observe  that  hierarchical  bases, 
while  being  very  natural  bases  for  the  p-version  stiffness  matrix,  are  veiy  unnatural  bases 
for  the  mass  matrix.  This  can  be  seen  by  noting  that  the  growth  in  the  condition  number 
is  exponential  in  p,  the  degree  of  the  polynomial  on  each  element 

We  often  need  to  solve  problems  using  the  mass  matrix,  in  particular,  time-dependent 
problems,  but  with  our  adaptive  codes  and  the  use  of  hierarchical  bases,  we  observe  that 
we  cannot  use  mass-lumping  to  remove  our  ill-conditioning. 

Using  methods  derived  from  those  for  the  stiffness  matrix,  we  shall  show  that 
numerically,  and  in  one  case  analytically,  we  can  control  this  ill-conditioning  at  little 
expense. 


J  DALIK 

An  approximation  of  the  two-dimensional  convection-diffusion 
problem  with  dominating  convection 

Let  us  consider  the  following  singularly  perturbed  boundary  value  problem: 

-  cAm  *  pVu  =  faj)  fa  Q  , 
a  -  u*  on  dQ  . 

Here  ft  is  a  bounded  domain  in  K2  and  e,  p  =  (p, ,  p2)  are  constant  functions  such 
that  0  <  e  <<  |p |  . 

We  present  a  numerical  method  relating  an  approximate  solution  uk  of  (1)  to  any 
triangulation  Tk  of  ft  with  the  discretization  step  h  greater  than  e .  This  method  can 
be  described  as  a  variant  both  of  the  finite-difference  method  and  of  the  Petrov-Galerkin 
method. 

We  formulate  the  following  two  results  of  theoretical  analysis:  First,  the  approximate 
solution  uk  does  not  oscillate.  Second,  the  discrete  L“-nomi  of  u  -  uk  on  certain 
subregions  5  £  ft  is  bounded  by  Ch2  +  Ke  .  This  local  error  estimate  is  proved  for 
strongly  regular  triangulations  only  and  is  uniform  whenever  no  layers  of  u  intersect 

S. 

By  means  of  graphical  illustrations  of  approximate  solutions  uk  of  the  problem  (1) 
and  of  its  nonstationary  equivalent,  we  indicate  that  the  breadths  of  layers  of  uk  depend 
on  the  Reynolds  number  |p| le  essentially  and  artificial  diffusion  has  no  influence  on 
them. 


L  DEMKOW1CZ  and  J  T  ODEN 

New  developments  on  application  of  hp-adaptive  BE/FE  methods  to  elastic  scattering 

Following  a  short  summary  of  the  ideas  and  results  described  in  [1,2,3]  two  new  aspects 
of  the  application  of  hp  approximations  to  the  solution  of  elastic  scattering  problems  are 
presented. 


28 


1.  Anisotropic  /(-refinement  techniques  and  the  related  issues  of  constrained 
approximation. 

2.  A  better  understanding  of  the  meaning  of  the  asymptotic  convergence  in  the 
context  of  the  problem  considered. 

A  number  of  typical  numerical  examples  conclude  the  presentation. 
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P  DRABEK 

Nonhomgeneous  degenerate  eigenvalue  problem 

Let  us  consider  the  nonhomogeneous  eigenvalue  problem 

-  <Mv(a(xji) I I*  *  Ab(jy()|ii|*',u  In  Q  , 

u  -  0  on  dO  .  (,b) 

We  assume  that  p  >  I  is  the  arbitrary  real  number  and  the  Carathlodory  functions 

i 

a(x,s),b(x,s)  satisfy  the  following  hypotheses.  Let  w  e  i^QX  —  €  L£l(Q)  a“d 

i  a(v)  i  cw(x)g(|*|)  , 

c  (2) 

0  <  Kv)  s  <*(*) 

for  a.e.  *  e  O  and  for  all  s  €  R,  where  c  >  0  is  a  constant,  a(x)  6  L"(fl)  and  g  :  R*  -* 
//,«)  is  a  nondecreasing  function. 

We  say  that  2  is  the  eigenvalue  and  u  €  WD'(w,C\)  is  the  corresponding  eigenfunction 
of  the  eigenvalue  problem  (1)  if  the  identity 


29 


f  a(xjt) | Vi» \[~*VuVWx  -  if  Hxrfluf  h&bc 

J  Q  JQ 


holds  for  any  <f>  £  W^fw.Cl)  (here  **  the  weighted  Sobolev  space). 

THEOREM  For  a  given  real  number  r  >  0  there  exists  the  least  eigenvalue  i,  >  0  and 
a  corresponding  eigenfunction  u,  a  0  (a.e.  in  A)  of  the  eigenvalue  problem  (1)  such  that 

”  r- 

Let  us  note  that  neither  global  results  for  nonlinear  eigenvalue  problems,  nor 
Ljusternik-Schnirelmann  theory,  can  be  used,  since  the  operators  are  not,  in  general, 
potential  operators. 

Let  us  note  also  that  the  principal  part  of  (1)  may  contain  the  degeneration  (or 
singularity )  and  it  may  depend  on  the  modulus  of  the  function  u  in  a  very  general  way 
(see  (2». 


M  G  EDWARDS 

A  flux  continuous  approximation  of  the  pressure  equation  for  h-adaptive  grids 

A  new  flux  continuous  discretisation  of  the  reservoir  simulation  pressure  equation  at  an 
h-adaptive  grid  interface  is  presented.  For  regular  grids  it  is  well  known  that  a  use  of 
a  harmonic  mean  of  the  permeabilities  to  define  the  discrete  pressure  equation 
coefficients  preserves  flux  continuity  and  is  important  in  the  case  of  discontinuous 
coefficients.  This  property  is  generalised  to  h-adaptive  grid  interfaces. 

The  standard  cell  centred  finite  volume  discretisation  of  the  pressure  equation  leads 
to  an  0(l/h)  pointwise  local  truncation  error  at  the  adaptive  grid  interface.  For  arbitrary 
grid  aspect  ratio  Ar,  and  interface  ratio  Gr,  it  has  been  shown  that  the  coefficient  of  this 
error  is  proportional  to  the  product  Ar‘Gr,  and  this  error  has  a  severe  impact  on  results 
for  large  Ar,  (1). 

A  correction  implemented  in  ( 1  ]  leads  to  a  non-symmetric  matrix,  conditional  diagonal 
dominance  and  increased  support  of  the  scheme  resulting  in  increased  complexity  in  the 
connectivity  of  the  pressure  matrix.  These  difficulties  are  handled  by  implementing  the 
correction  explicitly  following  (2]. 

In  this  paper  a  new  correction  is  presented  which  removes  the  local  leading  0(l/b) 
error.  The  new  correction  uses  the  same  support  as  the  standard  scheme  and  symmetry 
of  the  pressure  matrix  is  maintained.  The  correction  is  constructed  such  that  the  flux 
is  continuous  at  all  adaptive  grid  interfaces,  which  in  turn  results  in  an  unconditionally 
diagonally  dominant  pressure  matrix  for  isotropic  coefficients.  The  simplicity  of  the  new 
correction  enables  a  fully  implicit  and  consistent  implementation  of  the  scheme. 

Finally  the  benefits  of  the  new  correction  are  demonstrated  with  results  for  the  general 
case  of  a  dynamically  varying  h-adaptive  grid  which  tracks  a  moving  high  resolution 
shock  front 
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H  C  ELMAN 

Parallel  implementation  of  the  hp-version  of  the  finite  element  method 
on  a  shared-memory  computer 

We  study  the  costs  of  solving  two-dimensional  elliptic  partial  differential  equations 
discretized  by  the  hp-version  of  the  finite  element  method,  on  a  shared  memory  parallel 
computer.  The  computational  algorithm  consists  of  elimination  of  "high  order" 
unknowns  associated  with  element  interiors,  followed  by  preconditioned  conjugate 
gradient  solution  of  the  global  problem  on  element  boundaries,  using  the  submatrix 
associated  with  nodal  unknowns  as  a  preconditioner.  We  systematically  examine  costs 
in  CPU  time  of  the  individual  steps  of  the  computation,  including  construction  of  the 
local  stiffness  matrices.  Our  general  observations  are  that  costs  of  the  "naturally”  parallel 
computations  associated  with  local  elements  are  significantly  higher  than  any  global 
computations,  so  that  the  latter  do  not  represent  a  significant  bottleneck  to  parallel 
efficiency.  However,  in  order  to  avoid  inefficiency  caused  by  hierarchical  computer 
memories,  it  is  necessary  to  organize  the  local  operations  into  blocks  of  computations. 
We  show  how  to  do  this  effectively  using  the  Level  3  BIAS  linear  algebra  kernels. 


R  FALK 

Numerical  approximation  and  asymptotic  properties  of  some 
two  dimensional  plate  models 

The  finite  element  approximation  of  some  two  dimensional  plate  models  derived  from 
mixed  variational  principles  is  studied.  A  crucial  element  in  the  derivation  of  error 
estimates  is  to  first  understand  the  asymptotic  properties  of  the  models.  In  particular, 
we  consider  the  boundary  layer  behaviour  and  the  dependence  of  the  regularity  of  the 
solution  on  the  plate  thickness.  Differences  between  minimum  energy  models  and 
complementary  energy  models  are  highlighted. 


G  N  GATICA  and  G  C  HSIAO 

The  uncoupling  of  boundary  integral  and  finite  element  methods  for 
nonlinear  boundary  value  problems 


Recently,  the  usual  coupling  of  boundary  integral  and  finite  element  methods  has  been 
modified  to  yield  simpler  variational  formulations  for  linear  exterior  boundary  value 
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problems.  This  simplifed  procedure  is  based  on  choosing  the  artificial  coupling  boundary 
as  a  circle  or  a  sphere,  which  allows  us  to  invert  the  boundary  integral  operators  exactly. 
This  leads  to  the  so  called  uncoupling  technique,  by  means  of  which  the  resulting  weak 
formulations  reduce  to  almost  the  same  as  those  arising  from  Neumann  boundary  value 
problems  on  bounded  domains.  In  fact,  the  only  boundary  integral  operator  appearing 
in  the  formulations  is  the  one  provided  by  the  simple  layer  potential. 

The  purpose  of  this  paper  is  to  apply  the  uncoupling  procedure  to  study  the  weak 
solvability  of  certain  nonlinear  exterior  boundary  value  problems.  As  a  model  we 
consider  a  nonlinear  second  order  elliptic  equation  in  divergence  form  in  a  bounded 
inner  region,  which  becomes  the  Laplace  equation  in  the  corresponding  unbounded 
exterior  region.  We  provide  sufficient  conditions  for  the  nonlinear  coefficients  from 
which  existence,  uniqueness  and  approximation  results  are  established.  In  particular, 
nonlinear  equations  yielding  both  monotone  and  non-monotone  operators  are  analyzed. 
Also,  the  effect  of  numerical  integration,  and  the  polygonal  approximation  of  the 
coupling  boundary,  are  taken  into  account  to  perform  the  error  analysis  of  the  discrete 
scheme. 


L  GAVETE  CORVINOS,  C  MANZANO  DEL  MORAL  and  A  RUIZ  PEREA 
New  high  order  infinite  elements  with  different  decay  types 

A  class  of  problems  of  interest  in  mathematical  physics  and  engineering  is  characterized 
by  partial  differential  equations  on  unbounded  domains.  For  example,  the  irrotetional 
flow  of  an  incompressible  fluid  exterior  to  a  body  in  R’  is  described  by  Laplace’s 
equation  in  the  exterior  domain,  ft,  with  no  flow  through  the  body  boundary,  T,  and 
uniform  flow  at  infinity.  Similarly,  in  acoustics  and  electromagnetism,  Helmholtz’s 
equation  describes  the  exterior  scattering  from  a  body  in  R*.  In  this  case,  the  "boundary" 
condition  designates  the  decay  farfield  rate  of  the  solution  at  infinity.  Existence  and 
uniqueness  of  solutions  have  been  obtained  by  several  authors  using  various  techniques, 
but  the  treatment  of  unbounded  domains  has  presented  problems  for  the  finite  element 
analyst  Discretizations  are  typically  extended  to  large  distances  in  attempts  to  minimize 
inaccuracies.  Unfortunately,  coverage  of  this  extensive  domain  results  in  CPU  costs  and 
storage  penalities. 

The  focus  of  the  present  work  is  the  numerical  solution  of  partial  differential 
equations  by  using  infinite  elements  instead  of  finite  domain  truncation.  Infinite 
elements  are  becoming  increasingly  popular  as  an  economical  means  of  extending  the 
finite  element  method  to  deal  with  unbounded  domains.  Various  approaches  have  been 
adopted  in  the  extension  of  the  method,  and  these  techniques  tend  to  merge  into  each 
other.  In  this  paper,  new  high  order  mapped  infinite  elements  are  investigated. 

We  conclude  that  in  order  to  solve  efficiently  a  partial  differential  equation  with  an 
unbounded  exterior  domain,  it  is  necessary  to  consider  the  following: 

(a)  The  type  of  infinite  element  to  be  used  and  its  control. 

(b)  The  location  of  the  interface  between  finite  and  infinite  elements. 

(c)  Alternative  approaches  to  Finite  Element  Formulation  (only  in  some  cases). 


These  questions  can  be  solved  by  using  high  order  mapped  infinite  elements  with  the 
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appropriate  decay  type.  Also  we  consider  the  case  that  it  is  necessary  to  employ  a 
Petrov-Galerkin  method.  The  advantages  of  using  the  new  approach  proposed  in  this 
paper  are  shown  with  different  examples  of  solutions  of  partial  differential  equations  on 
unbounded  domains. 


R  GHANEM  and  J  E  AKIN 
A  stochastic  finite  element  analysis  method 

The  finite  element  method  is  <  n  appu*.d  to  partial  differential  equations  that  involve 
constitutive  relations  where  tL  sterial  constants  are  expected  to  have  random  vari¬ 
ations.  Certain  applications  are  c.  pec  ted  to  have  large  random  spatial  variations  in  the 
constitutive  parameters.  These  would  include  the  structural  or  electrical  response  of 
biological  materials  or  large-scale  environmental  problems,  such  as  flow  through  porous 
media.  This  paper  presents  the  extension  of  classical  finite  element  computational 
procedures  to  the  corresponding  complete  and  reliable  stochastic  finite  element  code. 

The  random  aspect  of  the  various  parameters  in  the  system  can  be  accounted  for 
through  the  introduction  of  an  additional  dimension  on  the  formulation  of  the  problem. 
Various  spectral  expansions  along  this  dimension  are  given  by  the  Karhunen-Loeve 
expansion  and  the  Polynomial  Chaos  expansion.  These  expansions  are  well  known  in 
applied  mathematics,  and  have  been  successfully  implemented  in  the  context  of  structural 
engineering  and  structural  mechanics.  These  expansions  are  incorporated  into  spatial 
discretizations  of  the  domain  in  accordance  with  standard  finite  element  procedures. 
The  Karhunen-Loeve  expansion  is  a  mean-square  convergent  optimal  expansion  of  a 
random  process  in  terms  of  a  denumerable  set  of  orthogonal  random  variables.  We 
provide  an  explicit  expression  for  the  solution  process  as  a  multidimensional  surface  in 
terms  of  orthogonal  polynomials,  the  coefficients  of  which  are  obtained  through  a 
Galerkin  projection  scheme  coupled  with  a  Polynomial  Chaos  expansion.  Such  an 
expression  is  superior,  both  in  accuracy  and  incomputationa!  efficiency,  to  currently  used 
first-  and  second-order  Taylor  expansions.  From  such  an  expansion,  the  probability 
distribution  function  of  the  solution  can  be  obtained. 

The  numerical  implementation  of  the  deterministic  projection  follows  the  standard 
guidelines  as  those  for  the  deterministic  finite  element  method.  The  treatment  in  this 
paper  will,  therefore,  be  confined  to  the  implementation  of  the  projection  in  the  space 
of  random  functions.  The  first  step  in  implementing  the  proposed  method  consists  of 
expanding  the  random  process  representing  the  spatial  randomness  of  the  system 
parameters  by  using  the  Karhunen-Loeve  expansion  of  the  process.  Then,  the 
computations  leading  to  the  expansion  consist  of  employing  an  integral  eigenvalue 
equation  that  can  usually  be  obtained  analytically.  The  next  step  in  the  process  consists 
of  implementing  the  Karbunen-Loeve  expansion  by  using  Polynomial  Chaos.  The  size 
of  the  resulting  system  of  equations  is  proportional  to  the  number  of  Polynomial  Chaoses 
used. 

The  technique  develops  a  fully  stochastic  finite  element  formulation  that  completely 
incorporates  the  stochastic  nature  of  the  problem.  The  outcome  is  an  accurate 
representation  of  the  solution  process  via  a  convergent  expansion  in  orthogonal 
stochastic  polynomials.  Unlike  first-  and  second-order  expansions,  this  representation 
affords  accurate  estimates  of  the  statistics  of  the  solution  beyond  the  second-order 
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moments.  These  statistics  can  be  either  in  the  form  of  higher-order  moments,  or  as 
convergent  approximations  to  the  probability  distribution  function  of  the  solution. 
Highly  accurate  approximations  of  the  probability  distribution  function  of  the  solution 
process  can  also  be  obtained. 


D  GIVOLI  and  J  B  KELLER 

l 

Non-reflecting  finite  elements 

Wave  problems  in  domains  extending  to  infinity  are  very  often  encountered  in  various 
fields  of  application  such  as  acoustics  and  structural  acoustics,  geophysics,  meteorology 
and  fluid  dynamics.  Standard  solution  procedures  of  such  problems  include  the 
boundary  element  method,  the  coupled  finite  element-boundary  element  method,  and 
the  use  of  infinite  elements.  Another  common  procedure  consists  of  the  following  three 
steps:  (a)  first  truncate  the  infinite  domain  by  introducing  an  artificial  boundary  B,  thus 
making  the  computational  domain  finite,  (b)  then  choose  a  certain  boundary  condition 
to  be  imposed  on  B,  (c)  finally  solve  the  problem  in  the  finite  computational  domain  by 
the  finite  element  method.  This  procedure  as  well  as  others  are  described  in  the  recent 
monograph  [1], 

It  is  well  known  (see  the  review  paper  [2])  that  a  poor  choice  of  the  artificial  boundary 
condition  on  B  may  give  rise  to  spurious  reflection  of  waves,  and  thus  to  large 
inaccuracies  in  the  finite  element  solution.  An  effective  boundary  condition  must  have 
the  property  that  outgoing  waves  hitting  B  are  transmitted  through  B  without  any 
reflection.  It  is  shown  in  (2]  that  a  robust  non-reflecting  boundary  condition,  effective 
in  the  general  case  and  not  only  in  some  specific  cases,  must  be  either  nonlocal  or  of 
sufficiently  high  order. 

Exact  nonlocal  boundary  conditions  were  devised  in  (3]  and  [4]  for  exterior  time- 
harmonic  acoustic  waves  and  elastic  waves.  The  artificial  boundary  B  was  chosen  to  be 
a  circle  in  two  dimensions  and  a  sphere  in  three  dimensions.  The  use  of  the  nonlocal 
boundary  conditions  leads  to  a  non-standard  finite  element  formulation.  The  infinite 
domain  outside  the  computational  domain  can  be  regarded  as  a  'super  finite  element* 
with  interacting  degrees  of  freedom  on  B.  The  finite  element  scheme  converges  to  the 
exact  solution  with  the  optimal  rate  of  convergence  (see  (2]  and  |3-7)). 

Local  high-order  non-reflecting  boundary  conditions  were  proposed,  e.g.  by  Engquist 
and  Majda  [8],  Bayliss  and  Turkel  {9]  and  Givoli  and  Keller  (4],  The  Bayliss-Turkel 
conditions  were  used  in  a  finite  element  scheme  in  [10].  The  use  of  boundary  conditions 
on  B  which  involve  high-order  spatial  derivatives  is  problematic  from  the  finite  element 
viewpoint,  since  standard  conforming  C°  finite  elements  cannot  be  used  in  that  case  in 
the  layer  of  elements  adjacent  to  the  boundary  B. 

la  the  present  paper  we  consider  both  nonlocal  and  local  non-reflecting  finite  elements 
(NRFEs),  and  discuss  the  differences  between  the  two  approaches  from  various 
perspectives.  We  devise  three  new  NRFEs: 

1.  A  nonlocal  NRFE,  used  with  an  elliptic  artificial  boundary.  This  NRFE  is 
particularly  useful  in  solving  scattering  problems  with  slender  obstacles. 

2.  A  nonlocal  NRFE  for  2D  problems  in  geophysics.  The  geometry  in  this  case  is 
semi-infinite,  and  the  construction  of  the  NRFE  turns  out  to  be  significantly  more 
difficult  than  in  the  case  of  exterior  problems. 


34 


3.  A  local  NRFE  with  high-order  continuity  on  the  boundary  B.  This  element 
possesses  C°  or  C2  continuity  on  one  of  its  sides,  and  C0  continuity  elsewhere. 
It  is  fully  compatible  with  standard  finite  element  architecture.  One  layer  of 
elements  of  this  type  is  used  near  the  artificial  boundary  0. 

We  will  discuss  the  computational  aspects  and  the  implementation  of  these  three 

elements,  and  will  demonstrate  their  performance  through  a  number  of  numerical 

examples. 
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G  A  GRAVVANIS  and  E  A  L1P1TAKIS 

Using  explicit  preconditioned  schemes  for  solving  initial/boundary  value  problems 
! 

A  new  class  of  hybrid  time  implicit-explicit  approximating  schemes  combined  with 
Approximate  Inverse  Finite  Element  Matrix  (AIFEM)  techniques  and  explicit 
prvconditioning  semi-direct  methods  for  the  numerical  solution  of  initial/boundaty  value 
problems  is  presented.  The  AIFEM  techniques  based  on  the  concept  of  adaptable  LDLT 
factorization  procedures  have  been  recently  introduced  for  computing  explicit 
pseudoinverses  of  large  sparse  symmetric  matrices  of  irregular  structure,  derived  from 
the  FE  discretization  of  Elliptic  and  Parabolic  Partial  Differential  Equations  without 
inverting  the  corresponding  decomposition  factors.  This  category  of  equations  represents 
a  large  class  of  commonly  occurring  problems  in  Mathematical  Physics  and  Engineering. 
The  AIFEM  techniques  originated  from  the  known  algorithmic  procedures  of  inverting 
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a  real  (n  x  n)  matrix  A,  which  can  be  decomposed  by  adaptive  approximate  factorization 
methods  (i.e.  A  =  L,DtL{,  where  D,  is  a  diagonal  matrix  and  L,  is  a  sparse  lower 
triangular  matrix  of  the  same  profile  as  A),  by  computing  explicitly  the  elements  of  the 
inverse  A'1,  without  inverting  the  decomposition  factors.  The  effectiveness  of  the  Explicit 
Preconditioned  iterative  methods  using  these  AIFEM  techniques  for  solving  numerically 
initial/boundaiy  value  problems  is  related  to  the  fact  that  the  exact  inverse  of  the  original 
sparse  coefficient  matrix  A  (although  full)  exhibits  a  similar  'fuzzy*  structure  to  A,  i.e. 
the  largest  elements  are  clustered  around  the  principal  diagonal  and  m-diagonal  where 
m  is  the  semi-bandwidth  of  A. 

It  should  be  noted  that  an  important  feature  of  the  AIFEM  algorithmic  techniques  is 
the  provision  of  both  explicit  preconditioning  direct  and  iterative  methods  for  solving 
partial  differential  equations  in  two  and  three  space  variables.  Additional  facilities  are 
provided  by  the  choice  of  'fill-in*  (factorization)  and  'retention*  (pseudoinverse) 
parameters  that  allow  the  best  method  for  a  given  problem  to  be  selected. 

Additionally  a  class  of  composite  Picard/Newton  methods  combined  with  Approximate 
Inverse  Finite  Element  Matrix  (AIFEM)  techniques  and  explicit  preconditioning  semi- 
direct  methods  for  the  numerical  solution  of  non-linear  elliptic  partial  differential 
equations  is  presented. 


D  V  GRIFFITHS 

Use  of  computer  algebra  systems  to  generate  element  stiffness  matrices 

Finite  element  matrices  such  as  those  for  stiffness  and  mass  are  usually  generated  using 
Gaussian  quadrature  because  it  leads  to  convenient  formulations  in  terms  of  local 
coordinates.  This  paper  describes  how  element  matrices  can  be  generated  in  dosed  form 
with  the  help  of  Computer  Algebra  Systems  (CAS)  and  bow  this  can  lead  to 
improvements  in  run-times.  The  CAS  approach  is  also  shown  to  be  able  to  generate 
matrices  which  would  normally  be  obtained  using  reduced  or  selective  reduced 
integration  (SRI). 

Consider  a  square  4-node  element  in  plane  strain,  fixed  on  three  sides,  with  a 
horizontal  force  P  applied  to  the  remaining  corner.  Computer  Algebra  can  be  used  to 
compute  the  horizontal  deflection  of  the  corner  for  various  element  stiffness  integration 
schemes.  The  expressions  show  why  locking  occurs  as  the  material  approaches 
incompressibility  with  full  integration. 
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C  GROSSMANN 

Newton’s  method  for  obstacle  problems 

The  obstacle  problem  can  be  considered  as  a  variational  problem  with  inequality 
constraints.  The  discretization  by  piecewise  linear  finite  elements  leads  to  large  scale 
optimizatibn  problems  with  a  special  structured  objective  functional  and  with  simple 
bounds  for  the  variable  as  constraints.  In  this  investigation  we  include  the  constraints 
into  an  auxiliaty  objective  functional  by  means  of  a  penalty  technique.  The 
unconstrained  problems  obtained  are  nonlinear  variational  equations.  These  can  be 
solved  by  Newton’s  method.  As  a  well  known  fact  the  auxiliary  problems  generated  in 
penalty  techniques  are  ill  conditioned  in  the  limit  Thus  an  important  question  is  to 
estimate  the  area  of  contraction  of  Newton's  method  in  dependence  of  the  penalty  and 
of  the  discretization  parameter. 

In  the  first  part  we  summarize  some  basic  facts  on  obstacle  problems,  their 
discretization  and  on  an  adapted  penalty  technique.  Later  we  derive  optimal  parameter 
selection  rules  to  adjust  the  error  caused  by  the  penalty  terms  to  the  same  magnitude 
as  the  discretization  error  of  the  finite  element  approximation.  Finally,  we  investigate 
the  convergence  behaviour  of  Newton's  method  applied  to  the  variational  equality 
obtained.  Here  the  area  of  contraction  of  Newton’s  method  is  estimated  in  dependence 
of  the  discretization  parameter  and  the  penalty  parameters  by  means  of  maximum 
principles  for  elliptic  problems.  Finally,  we  adjust  the  iterations  on  each  level  of 
discretization  and  investigate  the  method  on  families  of  grids. 
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1  F  GUARNACCIA  and  G  F  PINDER 

Animation  of  denser-than-water  immiscible  fluid  migration  in  the  subsurface 

Subsurface  migration  of  denser-than-water  immiscible  contaminants,  such  as  chlorinated 
hydrocarbons,  can  be  simulated  and  the  output  presented  via  animation.  The  four  non¬ 
linear  partial-differential  equations  that  describe  the  physics!  system  can  be 
approximated  using  orthogonal  collocation.  The  resulting  system  of  equations  are 
efficiently  solved  using  domain  decomposition  in  combination  with  Picard  iteration. 

Contaminant  saturation  and  concentration  are  effectively  represented  using  colour- 
based  computer  graphics  and  animation.  By  calculating  the  graphical  representation  of 
the  computed  values  of  these  parameters  for  each  time  step  of  the  simulation,  and 
converting  the  output  to  television  quality,  a  video  tape  of  the  dynamics  of  the  system 
is  generated.  Imbibition,  drainage  and  eventual  dissolution  of  the  contaminant  are 
represented  for  different  geohydrological  situations. 


J  GWINNER 

Boundary  element  methods  for  contact  problems  in  linear  elastostatics 

As  witnessed  by  the  monograph  of  Kikuchi  and  Oden,  [1],  contact  problems  in  nonlinear 
solid  mechanics  are  an  interesting  field  of  application  of  the  Finite  Element  Method. 
On  the  other  band,  contact  problems  show  much  more  concern  about  the  boundary 
behaviour,  namely  boundary  displacements  and  boundary  tractions  than  about  the  stress 
distribution  inside  the  body.  Therefore  it  is  natural  to  ask  for  the  application  of 
Boundary  Element  Methods,  at  least  for  linear  materials  within  an  infinitesimal  theory 
of  elastostatics.  There  is  recent  progress  towards  this  direction. 

In  this  contribution  we  address  the  most  interesting  case  from  the  view  of  mechanics 
and  the  most  delicate  case  from  the  view  of  mathematics  where  the  linear  elastic  body 
constrained  by  a  rigid  foundation  is  solely  subjected  to  applied  boundary  forces,  not 
necessarily  fixed  or  with  displacements  imposed  at  some  boundary  part  Using 
Somigliani's  identity  and  the  full  Calderon  projection  including  the  bypersingular  integral 
operator  a  boundary  integral  variational  inequality  [3]  can  be  derived  as  a  mixed 
variational  formulation  for  the  boundary  displacements  and  boundary  tractions.  The 
resulting  bilinear  form  is  not  symmetric,  but  nonnegative  and  semicoercive  in  the  sense 
that  it  satisfies  a  Garding  inequality  [3]. 

Use  of  periodic  spline  functions  as  trial  and  test  functions  on  the  boundary  leads  to 
Galerkin  boundary  element  methods.  Here  we  consider  not  only  piecewise  linear 
approximations  for  the  boundary  displacements  and  piecewise  constants  for  the  boundary 
tractions,  but  also  approximations  of  higher  order  that  lead  to  a  nonconforming 
approximation  of  the  unilateral  constraint  of  the  contact  problem.  Based  on  a  recent 
descretization  theory  for  semicoercive  variational  inequalities  [2]  convergence  of  those 
boundary  element  methods  can  be  established  with  respect  to  the  energy  norm  [3]. 

Further  we  report  ou  numerical  calculations  [4]  for  the  well-known  Hertz  contact 
problem.  By  a  Schur  complement  technique  we  transform  the  discrete  variational 
inequality  into  a  linear  complementarity  problem  that  is  formulated  only  in  terms  of 
displacements.  As  an  iterative  solver  the  projected  SOR  algorithm  is  investigated. 
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Moreover  in  virtue  of  a  recent  extension  of  tbe  discretization  theory  to  variational 
inequalities  of  the  second  kind  [5]  we  can  also  treat  contact  problems  with  given  friction. 
To  discretize  the  nonlinear  friction  functional  various  quadrature  rules  can  be  employed 
that  can  be  analysed  in  a  similar  way  as  for  the  finite  element  method  (6j.  Finally  we 
discuss  bow  more  realistic  friction  models  can  be  incorporated  in  the  boundary  integral 
approach. 
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S  HAMMARLING 

LA  PACK:  High  performance  software  for  linear  algebra 

NAG  Ltd  and  US  Scientists  have  collaborated  in  tbe  development  of  LAPACK,  a  high 
performance  linear  algebra  package,  already  proven  for  shared  memory  parallel 
machines,  now  being  extended  to  distributed  memory  machines. 


I  HARArtl  and  T  J  R  HUGHES 

Finite  element  methods  for  general  second-order  elliptic  partial  differential  equations:  Model 
problems 

Finite  element  methods  for  second-order  elliptic  partial  differential  equations  were 
initially  applied  to  problems  which  contain  only  derivatives  of  the  highest  order.  In  such 
cases  finite  elements  based  on  the  Galerkin  method  possess  several  advantageous 
attributes,  with  guaranteed  performance  on  configurations  of  practical  interest.  The 
addition  of  other  terms  such  as  first-order  derivatives  in  advective-diffusive  problems  and 
undifferentiated  terms  in  the  Helmholtz  equation  has  potentially  destabilizing  effects. 
The  concept  of  Galerkin/least  squares  (GLS)  obtained  by  appending  terms  in  least 
squares  form  to  the  standard  Galerkin  formulation  enables  the  design  of  stable  methods 
for  such  applications,  retaining  the  higher-order  accuracy  of  tbe  Galerkin  method  for 
problems  governed  by  the  Laplacian,  at  relatively  low  added  computational  cost. 

General  second-order  equations,  which  are  studied  in  this  work,  have  exact  solutions 
which  may  differ  widely  in  nature  depending  on  values  of  the  coefficients,  with  rapidly 
varying  values  in  domain  interiors  as  well  as  close  to  boundaries.  An  analysis  of 
numerical  solutions  obtained  by  tbe  Galerkin  method  employing  linear  finite  elements 
indicates  that  relatively  fine  meshes  may  be  required  to  retain  an  acceptable  degree  of 
accuracy  for  certain  ranges  of  values  of  the  physical  coefficients.  Numerous  criteria  for 
improving  accuracy  by  employing  GLS  technology  are  examined.  Of  these,  several  are 
found  to  provide  the  lowest  error  in  nodal  amplification  factors,  each  in  a  certain  range 
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of  ratios  of  physical  coefficients.  Guidelines  for  required  mesh  resolutions  are  presented 
for  Galerkin  and  GLS,  showing  that  in  large  parts  of  the  parameter  space  substantial 
savings  may  be  obtained  by  employing  GLS.  Numerical  computations  bear  out  these 
observations. 


B  HEINRICH 

The  Fourier-finite  element  method  for  elliptic  problems  with  axisymmetric  edges 

The  Fourier-finite  element  method  is  a  semianatytic  technique  for  the  approximate 
solution  of  three-dimensional  elliptic  problems  in  axisymmetric  domains  with  non- 
axisymmetric  data.  The  method  employs  truncated  Fourier  series  with  coefficients  which 
are  the  solutions  of  a  finite  number  of  two-dimensional  elliptic  problems  posed  on  the 
meridian  plane  of  the  domain  and  approximated  by  the  finite  clement  method.  This 
approach  is  often  used  by  engineers  for  thermal  and  stress  analysis  of  axisymmetric 
bodies  and  can  be  implemented  on  a  parallel  computer. 

In  this  paper  we  present  results  on  the  mathematical  justification  of  the  Fourier-finite 
element  approach  for  solving  problems  in  3D,  especially  for  Poisson’s  equation,  the 
description  of  the  properties  of  the  approximation,  local  mesh  refinement  in  the 
meridian  plane  for  an  appropriate  treatment  of  singularities  near  re-entrant  edges  and 
error  estimates  in  the  Hl-  and  Lg-norms,  with  respect  to  the  discretization  parameters 
N  and  h  (Fourier  and  finite  element  approximation)  for  solutions  belonging  to  Sobolev 
spaces.  Further,  algorithmic  aspects  and  the  computer  implementation  of  the  method 
are  discussed. 


J  P  HENNART 

Nodal  finite  element  methods  for  partial  differentia!  equations 

Modern  nodal  methods  were  developed  in  the  1970's  in  numerical  reactor  calculation. 
Roughly  speaking,  nodal  methods  are  fast  and  accurate  methods  which  try  to  combine 
the  attractive  features  of  the  finite  difference  method  (FDM)  and  the  finite  element 
method  (FEM).  From  the  FEM,  they  borrow  a  piecewise  continuous,  usually  (but  not 
always)  polynomial,  behaviour  over  a  given  coarse  mesh.  With  the  FDM,  they  have  in 
common  the  fact  that  the  final  algebraic  systems  are  usually  quite  sparse  and  well 
structured,  at  least  over  domains  which  are  not  too  irregular  such  as  unions  of  rectangles. 
Nodal  methods  can  thus  be  viewed  « » fast  solvers,  to  which  techniques  of  vectorization 
and  (or)  parallelism  are  applicable.  They  are  especially  suited  to  all  physical  situations 
which  have  been  traditionally  modelled  by  finite  differences,  and  for  which  the  domain 
as  well  as  the  coefficients  of  the  equation(s)  are  not  known  with  sufficient  accuracy.  It 
is  then  tempting  to  use  a  rectangular  grid  to  discretize  the  domain  and  to  assume  that 
over  each  cell  of  the  grid  the  physical  parameters  are  not  known  with  great  accuracy,  so 
that  a  constant  or  mean  value  only  is  available.  This  situation  is  prevailing  for  instance 
in  groundwater  hydrology,  oil  reservoir  simulation,  and  air  pollution  modelling  problems. 

In  this  paper,  we  show  that  most  nodal  schemes  can  be  viewed  as  finite  element 
schemes  provided  the  FEM  is  viewed  in  a  fairly  liberal  way  with  moments  (and  not  point 
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values)  as  principal  unknowns,  with  piecewise  polynomials  (or  quasi-polynomials)  as 
basic  functions  and  with  a  systematic  use  of  variational  crimes  like  numerical  quadrature 
and  nonconformity. 

These  ideas  were  applied  to  the  three  basic  types  of  PDE’s,  elliptic,  parabolic,  and 
hyperbolic,  with  concrete  applications  (and  numerical  results)  in  multidimensional 
diffusion,  space-time  kinetics  and  2D  transport  of  neutrons  in  nuclear  reactors,  as  well 
as  for  benchmark  problems  of  the  types  mentioned  above,  including  the  Poisson  equation 
and  the  heat  diffusion  equation. 


I  HERRERA 

Localized  Adjoint  Method:  An  overview 

The  Localized  Adjoint  Method  (LAM)  is  a  new  and  promising  methodology  for 
discretizing  partial  differential  equations,  which  is  based  on  Herrera’s  Algebraic  Theory 
of  Boundary  Value  Problems.  Applications  have  successively  been  made  to  ordinary 
differential  equations,  for  which  highly  accurate  and  efficient  algorithms  were  developed, 
multidimensional  steady  state  problems  and  optimal  spatial  methods  for  advection- 
diffusion  equations.  More  recently,  generalizations  of  Characteristic  Methods  known  as 
Eulerian-Langragian  Localized  Adjoint  Method  (ELLAM),  were  developed  and  many 
specific  applications  have  already  been  made.  Related  work  and  additional  applications 
are  underway. 

In  the  construction  of  approximate  solutions  there  are  two  processes,  equally 
important  but  different.  They  are: 

i)  Gathering  informaiton  about  the  sought  solution;  and 

ii)  Interpolating  or,  more  generally,  processing  such  information. 

These  two  processes  are  distinct,  although  in  many  numerical  methods  they  are  not 
differentiated  clearly.  The  information  about  the  exact  solution  that  is  gathered,  is 
determined  mainly  by  the  weighting  functions,  while  the  manner  in  which  it  is 
interpolated  depends  on  the  base  functions  chosen.  Examples  have  been  given  for  which 
these  processes  are  not  only  independent  but,  they  do  not  need  to  be  carried  out 
simultaneously. 

The  questions  posed  by  the  above  comments  are  quite  complex  and  to  explore  them 
in  all  their  generality  is  very  difficult  Localized  Adjoint  Method  is  a  methodology  we 
have  proposed  to  carry  out  such  analysis  and  develop  new  numerical  procedures  using 
the  insight  so  gained. 

In  this  paper  Localized  Adjoint  Method  is  briefly  explained  and  ELLAM  procedures, 
which  have  been  quite  successful  for  treating  advection  dominated  transport,  are 
discussed. 


Z  KESTftANEK 

Finite  elements  for  extended  variational  formulation  in  nonlinear  problems 

Nonlinear  problems  such  as  the  elastohydrodynamic  lubrication  problem,  involving 
nonconservative  mechanical  systems  are  of  considerable  importance  in  the  field  of  hydro- 
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and  aero-elasticity  and  space  mechanics.  These  problems  involve  nonsymmetric 
operators,  and  therefore  variational  formulations  in  a  strict  sense  are  not  available.  The 
variational  formulation  presented  [1],  (2]  is  called  extended  in  the  following  sense: 
Given  a  problem  N(q)  =  0,  with  a  nonlinear  operator  N  :  D(N)  C  U  -*  R(N)  C  V  = 
U\  find  a  functional  F,  if  any,  whose  critical  points  are  solutions.to  the  problem  and  yice 
vgrsa.  This  implies  that  for  a  given  operator  .N  an  operator  N  exists  such  that  5F  = 
(N(u),5uN)  and  the  problems  N(u)  =  0V  and  N(u)  *  0,  have  the  same  solution.  This 
statement  requires  only  that  the  critical  points  should  coincide  with  the  solutions,  without 
posing  the  additional  requirement  that  N  must  be  the  gradient  of  the  functional.  The 
problem  can  be  rewritten  as  N(u)  =  N;*(u;KN(u))  =  0,  provided  that  D(K)  C  R(N), 
R(K)  C  Q(N,*),  where  K  is  a  linear,  invertible  and  symmetrical  operator..  The.modified 
operator  N  satisfies  the  symmetry  condition  foe.  the  Gateaux  derivative,  N;  =N;*  and  it 
follows  that  the{e  exists  a  potential  functional  F[u]  =  1/2  (N(u),  KN(u))  whose  gradient 
is  the  operator  N.  It  is  worth  noting  that  the  operator  K  is  not  uniquely  determined  and 
the  practical  construction  of  K  for  the  general  case  presents  some  difficulties.  In 
numerical  applications,  however,  K  can  easily  be  defined.  The  functional  vanishes  when 
the  solution  is  reached.  Moreover  if  K  is  positive  definite,  then  F[u]  is  minimum  at  the 
critical  point  Functional  F  can  be  defined  over  finite-dimensional  subspaces  S.  The 
kernel  k(s,t)  of  the  linear  operator  K  in  integral  form  is  naturally  defined  over  S  x  S  and 
may  be  expressed  as  k(s,t)  =  (L(s)}T[A]{L(t)},  where  [A]  is  a  symmetric  constant  matrix 
and  {L}  is  a  vector  of  the  shape  functions.  Here  the  choice  of  the  kernel  in  a  finite 
subspace,  generated  by  means  of  functions  with  local  compact  support,  makes  it  possible 
to  provide  finite  element  models  in  a  straightforward  way.  The  main  novelty  in  this 
standard  procedure  is  a  complication  arising  from  the  intcgro-differential  nature  of  the 
basic  equation.  To  find  the  minimum  of  the  functional  F  we  can  use  the  Newton-like 
type  algorithm.  This  method  for  the  solution  of  the  inverse  problem  has  been 
successfully  developed  for  special  numerical  applications. 
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AQM  KHALIQ 

A  parallel  semidiscrete  Galerkin  splitting  method  for  second-order  hyperbolic  equations 

Second-order  hyperbolic  partial  differential  equations  (PDEs)  in  multi-space  dimensions 
occur  frequently  in  viscoelastic  theory  [3],  acoustics  [4],  and  in  many  geoscience 
applications,  for  example  seismology  (1,2, 5).  Computationaltechniquesareproposedby 
several  authors,  see  for  example  [2, 3, 6, 7]  and  references  therein,  for  the  numerical 
solution  of  the  wave  equation  by  finite-element  methods.  Their  investigations  have  been 
restricted  to  mainly,  two-dimensional  wave  equations  and  implementation  on  aerial 
machines.  Recently  an  explicit  finite-difference  scheme  has  been  propposed  by  Mufti, 
[5),  for  the  numerical  solution  of  large-scale  three-dimensional  seismic  models.  Since 
analytical  methods  are  also  not  effective  because  they  are  restricted  to  simple  geometries 
with  homogeneous  structures  like  the  finite-difference  methods,  there  is  an  ever  growing 
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need  to  develop  efficient  and  stable  finite-element  methods  which  may  be  implemented 
effectively  on  modern  parallel  computers. 

A  new  parallel  method  for  solving  hyperbolic  PDEs  in  many  space  dimensions  is 
presented  by  applying  the  finite-element  technique,  and  is  unconditionally  stable.  The 
proposed  method  is  based  on  the  semidiscrete  Galerkin  approach  utilizing  Strang-type 
splitting  techniques  in  which  a  multi-dimensional  problem  is  reduced  to  a  series  of  one¬ 
dimensional  problems.  This  approach  avoids  the  difficulty  of  constructing  a  finite- 
element  space  for  multi-dimensional  problems.  A  large  scale  problem  requiring  large 
storage  and  CPU  times  is  thus  solved  using  only  a  one-dimensional  finite-element 
method  on  MIMD  parallel  machine.  This  approach  is  seen  to  be  highly  competitive  to 
those  presented  in  [7].  Second  order  convergence  and  error  estimates  of  the  scheme  are 
given  using  piecewise  linear  basis  functions.  However,  it  is  anticipated  that  cubic  B~ 
splines  may  also  be  used  as  basis  functions  and  their  convergence  may  be  proved  by 
semi-group  theory.  Numerical  results  are  demonstrated  on  problems  from  the  literature. 
The  proposed  parallel  algorithm  for  the  numerical  solution  of  the  wave  equation  will 
provide  geophysicists  with  a  very  useful  and  efficient  tool  for  predicting  and 
understanding  seismic  wave  propagation.  It  will  also  be  useful  to  engineers  because  of 
its  efficiency  for  solving  large  scale  problems  in  acoustics  and  in  numerical  approximation 
to  the  solution  of  the  linear  elastodynamic  equations. 
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B  KOVACS  and  F  J  SZAB 6 
On  the  stability  of  layered  circular  rings 

The  purpose  of  this  paper  is  to  find  the  critical  external  pressure  of  layered  circular  ring 
bucklings  in  their  planes.  The  composite  ring  subjects  to  normal  gas  pressure  and 
constant-directional  pressure.  The  stability  of  an  arch  is  analysed  by  formulating  the 
governing  differential  equations  and  the  associated  boundaty  conditions  using  a 
variational  principle.  The  cross  sections  in  individual  layers  remain  planes  after 
deformation. 

In  the  case  of  free  rings,  the  critical  external  pressures  are  investigated,  using 
theoretical  method,  the  COSMOS/M  and  the  SYSTUS  finite  element  software  packages. 
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As  a  theoretical  investigation  the  eigenproblem  is  solved  by  a  successive  approximation 
based  on  the  Schwarz-method  and  by  the  Rite-method,  reducing  the  problem  to  a 
matrix-eigenvalue  problem. 

For  the  finite  element  solution  the  model  of  the  ring  is  built  up  using  four-nodes 
composite  shell  elements  and  the  numerical  results  are  obtained  by  the  Inverse  Power 
iteration  technique  in  case  of  COSMOS/M  and  by  the  Lanczos  algorithm  in  SYSTUS. 

The  analysis  of  the  numerical  results  shows  that  in  the  case  of  free  rings  the  constant- 
direcitona!  pressure  gives  higher  critical  values  than  the  hydrostatic  pressure  and  one  can 
find  a  good  agreement  between  the  theoretical  and  finite  element  solutions. 


B  KOVACS  and  F  J  SZAB6 

Free  vibrations  of  layered  circular  ring  segments 

The  free  vibration  eigenvalue  problem  of  circular  arc  shaped  three  layered  sandwich  ring 
segments  was  investigated  theoretically  and  numerically  in  the  case  of  linear  clastic  layer- 
materials.  The  comparison  of  the  different  solution  methods  makes  it  possible  to  study 
the  plane  bending  free  vibration  behaviour  of  the  curved  layered  beams  more  accurately. 

The  free  vibration  eigenvalue  problem  of  the  ring  segments  was  solved  theoretically 
by  using  the  Schwarz-method,  which  is  a  fast  and  flexible,  easy-to-program  algorithm. 
Ihis  algorithm  can  be  easily  completed  by  an  optimization  method  and  can  be  applied 
for  some  special  optimization  problems,  which  cannot  be  handled  using  the  built-in 
optimization  of  the  known  finite  element  programs  (multiobjective  optimization,  discrete 
variables,  complicated  constraints,  etc.). 

The  results  of  the  theoretical  solution  were  compared  with  those  obtained  by  a  finite 
element  analysis  made  in  the  SYSTUS  finite  element  program  which  runs  on  a  VAX 
3100  workstation  under  VMS.  For  the  eigenvalue  solution  the  LANCZOS  method  was 
used.  Two  different  finite  element  models  were  built  up: 

The  first  model  contains  eight-nodes  hexahedron  brick  elements  and  this  makes 
it  possible  to  apply  several  elements  along  the  thickness,  which  increases  the  accuracy 
of  the  model. 

The  second  model  consists  of  four-nodes  composite  shell  elements,  which  results 
in  a  very  easy-to-solve  finite  element  model  which  needs  much  less  storing  and 
computation  capacity  than  the  previous  one. 

Numerical  results  of  the  theoretical  and  finite  element  analysts  were  made  for  the 
circular  arc  of  radius  varied  in  the  ranges  between  100mm  to  200mm,  with  the  same 
thickness  and  width  data,  in  the  case  of  two  different  core  materials. 

(Although  both  of  the  first  four  bending  eigenvalues  shows  a  very  close  agreement 
between  the  theoretical  and  finite  element  results,  which  is  better  in  case  of  long  (length 
of  150mm  or  more)  beams.  In  the  case  of  very  short  (length  of  50mm  or  less)  beams 
there  may  be  a  considerable  difference  between  the  results,  so  for  these  structures  it 
could  be  necessary  to  investigate  the  problem  experimentally,  too.  Comparing  the  two 
different  finite  element  models,  one  can  say  that  the  model  containing  hexahedron 
elements  needs  more  computing  and  storing  capacity  but  gives  more  accurate  results 
than  composite  shell  elements.  Because  of  the  possibly  long  calculation  time  in  case  of 
complicated,  long  solution  procedures  (nonlinear or  transient  dynamic  analysis  containing 
many  time  cards)  the  composite  shell  elements  are  more  suitable. 
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M  KftteEK 


Finite  element  approximation  of  a  nonlinear  anisotropic  heat  conduction  problem 

We  prove  the  existence  and  uniqueness  of  a  solution  of  the  following  second  order 
quasilinear  elliptic  problem  of  a  divergence  form  with  nonlinear  mixed  boundary 
conditions: 


divC4(x^<)grad  u)  -  F(xM) 

in  SI , 

U  m  U 

on  r, , 

nrA(sM)gnd  u  *  G(sj »)  -  0 

on  r,. 

where  ftCR'.de  {1,2,...},  is  a  bounded  domain  with  a  Lipschitz  boundary  afl,  n  is  the 
outward  unit  normal  to  the  boundary,  T,,  P2  are  relatively  open  and  disjoint  subsets  of 

afl  and  ?,  vj  Fj  -  30,  f  €  wftO)  a ad  A  *  («„)^. ,  «s  *  uniformly  positive  definite 

matrix  function.  We  further  assume  that  F(x,.)  and  G(r,.)  are  nondecreasing  functions 
for  almost  every  r£fl  and  s  €  T2  (in  particular,  they  can  be  independent  of  u).  We 
show  that  the  problem  is  nonpotential  in  general  and  also  that  the  theory  of  monotone 
operators  cannot  be  applied.  Sufficient  conditions  from  [1]  which  guarantee  the 
uniqueness  of  the  classical  and  weak  solutions  are  introduced.  Note  that  there  exist 
examples  of  non-unique  solutions  if  the  equation  is  not  in  divergence  form  (see  e  g.  [2]). 
We  also  gjve  a  finite  element  approximation  of  the  problem  and  present  some  results 
concerning  the  existence  and  uniqueness  of  the  Galerkin  solution.  We  then  introduce 
two  convergence  theorems  of  (1).  A  practical  example  of  computing  a  temperature  field 
in  the  magnetic  cores  of  transformers  is  shown. 
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M  S  KUCZMA  and  J  R  WHITEMAN 
A  variational  inequality  formulation  for  flow  theory  plasticity 

We  consider  the  quasi-static  deformation  process  of  an  elasto-plastic  body,  under  the 
assumption  of  small  strains.  The  elasto-plastic  behaviour  of  the  material  is  assumed  to 
be  governed  by  the  von  Mises  yield  function  with  linear  strain  hardening  (softening). 
The  yield  function  is  written  in  terms  of  strains,  and  we  define  the  incremental  problem 
as  the  system  that  consists  of  (a)  a  variational  equation  being  the  equilibrium  condition 


for  the  body  and  (b)  a  variational  inequality  expressing  the  unilateral  character  of  the 
plastic  strain-rate  multiplier.  An  iterative  method  for  solving  this  non-linear  system  is 
proposed,  and  results  of  some  numerical  experiments  based  on  the  FEM-approximation 
are  provided. 


G  KUHN 

A  BE-formulation  for  finite  strain  elastoplasticity 

An  extension  of  the  boundary-element-method  to  finite  strains  and  to  constitutive 
relations  based  on  the  concept  of  an  intermediate  configuration  is  presented.  Within  this 
concept  material  behaviour  is  described  in  terms  of  a  hyperelastic  relation  relative  to  the 
intermediate  configuration  together  with  evolution  equations  for  this  configuration  and 
further  internal  variables.  In  the  context  of  numerical  solution  techniques  this  enables 
an  operator  split  leading  to  a  global  hyperelastic  problem  and  local  time  integration 
schemes  for  the  internal  variables.  The  proposed  BEM-approach  allows  for  arbitrary 
hyperelastic  relations  relative  to  the  intermediate  configuration  and  leads  to  a  nonlinear 
set  of  equations  with  displacement  gradients  as  basic  unknowns.  As  a  Total-Lagrange 
scheme  is  adopted  the  expensive  computation  of  the  system  matrices  has  to  be  done  only 
once  in  the  initial  configuration.  Nevertheless  the  nonlinear  set  of  equations  can  be 
consistently  linearised  with  respect  to  the  local  integration  schemes.  Constitutive 
equations,  integration  schemes  and  linearisation  techniques  used  in  FEM  can  be  directly 
applied.  The  BEM-formulation  is  presented  for  associated  rate  independent 
elastoplasticity  using  a  yield  condition  in  strain  space.  Some  numerical  results  will  be 
shown. 


G  KUNERT 

Adaptive  mesh  refinement  in  the  finite  element  method 

The  rate  of  convergence  of  the  finite  element  method  can  be  improved  by  applying 
adaptive  mesh  refinement.  This  has  been  investigated  for  different  kinds  erf  error 
indicators,  and  numerical  experiments  were  performed  for  the  Laplace  operator. 


R  LEVKOVITZ  and  G  MITRA 

The  factorization  of  unstructured  sparse  symmetric  positive  definite  matrices 
on  distributed  memory  MIMD  parallel  computers 

The  efficient  solution  of  unstructured  large  sparse  symmetric  positive  definite  equations 
on  a  distributed  memory  computer  using  direct  methods  is  an  important  and  challenging 
area  of  parallel  computing  research.  We  concentrate  on  the  design  and  implementation 
of  a  parallel  Cholesky  factorization  algorithm  and  focus  the  discussion  on  the  distribution 
of  data,  efficient  distribution  of  the  computation  and  the  aggregation  of  the 
communication.  We  show  that  the  accumulation  of  non-zeros  in  certain  areas  of  the 
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Cholesky  factor  tends  itself  to  a  hybrid  sparse-dense  computation  paradigm  and 
demonstrate  how  to  extend  this  approach  to  a  distributed  memory  parallel  MIMD 
computer. 


J  MACKENZIE,  T  SONAR  and  E  SULI 
Adaptive  finite  volume  methods  for  hyperbolic  problems 

We  develop  an  a  posteriori  error  analysis  of  finite  volume  methods.  Local  a  posteriori 
error  estimates  are  derived  in  various  norms  which  provide  local  bounds  for  the  global 
error  in  terms  of  a  computable  residual.  For  hyperbolic  problems  with  smooth  solutions 
the  error  and  the  residual  are  measured  in  local  L2-norms;  however,  when  the  solution 
is  discontinuous,  as  is  typically  the  case  for  nonlinear  hyperbolic  conservation  laws,  we 
use  a  local  negative  Sobolev  norm,  instead. 

We  apply  these  residual-based  a  posteriori  error  estimates  to  two  finite  volume 
methods.  The  first  one  is  the  cell  vertex  finite  volume  method  on  structured 
quadrilateral  meshes  [1].  The  numerical  investigation  of  the  performance  of  the  error 
indicator  for  a  model  problem  of  two-dimensional  advection  demonstrates  its  reliability 
and  efficiency  [2].  The  error  indicator  is  then  applied  to  the  transonic  Euler  equations 
where  it  is  contrasted  with  classical  techniques  based  on  calculating  gradients  of  flow 
variables. 

The  performance  of  one  of  the  residual-based  error  estimators  for  a  two-dimensional 
linear  advection  problem  solved  by  the  cell  vertex  method  is  shown.  Very  good 
agreement  can  be  found  between  the  error  estimate  and  the  local  error. 

The  secbnd  adaptive  finite  volume  method  is  constructed  on  a  general  triangulation, 
and  it  is  based  on  the  Riemann  solver  of  Osher  and  Solomon  [3].  A  MUSCL  recovery 
technique  is  employed  to  increase  the  accuracy  through  baiycentric  subdivisions.  The 
adaptive  algorithm  uses  point  insertion  techniques.  Again,  the  finite  element  residual 
serves  as  an  error  indicator  to  control  the  adaption  process.  A  triangle  K  of  the 

triangulation  is  refined  if  >  TOL  or,  alternatively,  if  |r*|>  tw  >  TOL. 

Our  numerical  results  indicate  the  potential  of  residual-based  a  posteriori  error 
estimates  for  finite  volume  approximations  of  compressible  flow  problems  and  highlight 
their  superiority  when  compared  with  classical  methods  of  adaptive  error  control. 
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H  MAISCH,  H  KARL  and  G  LEHNER 

The  application  of  the  Finite  Element  Method  to  the  solution  of  nonlinear  partial  differential 
equations  and  a  comparison  to  the  Finite  Difference  Method. 

Solitons  are  special  solutions  of  certain  nonlinear  partial  differential  equations  with 
increasing  importance  in  plasma  physics,  hydrodynamics,  electrical  engineering  and  other 
fields.  In  this  paper  we  have  studied  the  applicability  of  the  finite  element  method  to 
solve  soliton  equations.  As  the  most  important  and  famous  example  we  have  solved  the 
Korteweg-de  Vries  equation  (KdV)  using  Hermite  elements.  Here  we  have  studied  the 
interaction  of  solitons  among  themselves  and  also  the  generation  of  solitons  from 
arbitrary  initial  conditions.  The  results  agree  very  well  with  theoretical  predictions.  In 
the  case  of  soliton  generation  there  are  often  no  analytical  results  available,  therefore 
we  have  checked  the  three  lowest  order  conservations  laws  While  for  the  Korteweg-de 
Vries  equation  the  continuity  of  the  first  derivative  is  absolutely  necessary,  for  most 
other  equations  continuity  conditions  are  less  stringent  and  Lagrange  type  elements  are 
sufficient.  Some  numerical  results  for  the  regularized  long  wave  equation  (RLW),  the 
Sine-Gordon  equation  (SG),  nonlinear  Schrddinger  equation  (NLS)  and  the  Zakharov 
equations  are  also  presented  briefly.  The  numerical  solution  of  two  dimensional  soliton 
equations  is  now  also  possible  due  to  increasing  computer  power.  We  have  chosen  the 
two  dimensional  Sine-Gordon  equation  as  an  example.  For  comparison  we  solved  the 
Korteweg-de  Vries  equation  and  the  two-dimensional  Sine-Gordon  equation  also  by 
using  finite  differences,  especially  the  method  of  lines.  It  turned  out,  that  the  more 
flexible  finite  element  method  definitively  performed  better  in  both  cases. 


J  MANDEL 

Balancing  domain  decomposition  preconditioners  and  discontinuous  coefficients 

This  talk  will  present  the  formulation  and  experience  with  several  preconditioners  based 
on  solving,  in  each  iteration,  a  local  problem  on  each  subdomain  coupled  with  a  coarse 
level  problem.  This  provides  for  an  efficient  global  propagation  of  the  error  and 
guarantees  that  the  possibly  singular  local  problems,  solved  in  every  iteration,  are 
consistent. 

The  preconditioner  is  based  on  solving  in  each  iteration  possibly  singular  problems 
based  on  the  local  stiffness  matrices  of  the  subdomains.  Is  it  proved  that  the  condition 
number  remains  bounded  for  arbitrary  size  jumps  of  coefficients  between  subdomains 
in  any  number  of  dimensions,  is  bounded  independently  of  the  number  of  subdomains, 
and  grows  only  as  the  logarithm  of  subdomain  size.  Computational  experiments  confirm 
the  theory  and  show  that  the  method  is  remarkably  robust  and  performs  very  well  for 
general  strongly  dicontinuous  coefficients  as  well  as  unstructured  subdomains.  Possible 
modifications  for  shells  and  plates  will  be  also  discussed. 


J  MANDEL 


Fast  iterative  solvers  for  p-version  solids,  plates,  and  shells 

We  present  the  formulation  of  and  experience  with  several  methods  for  the  iterative 
solution  of  large-scale  systems  of  algebraic  linear  equations  arising  from  the  p-version 
finite  element  method  in  three  dimensions.  The  basic  idea  of  the  method  is  to  use 
domain  decomposition  approaches  with  each  element  considered  to  be  a  subdomain. 
The  resulting  iterative  methods  have  been  shown  to  be  superior  to  state-of-the-art  direct 
solvers  for  real-world  problems  with  severely  distorted  solid  elements,  such  as  engine 
crankshaft,  complete  aircraft  fuselage,  and  stiffened  shell  panel  with  a  window  rim.  The 
data  have  been  generated  by  the  packages  MSC/PROBE  and  STRIPE. 

The  method  is  essentially  block  preconditioning,  with  one  block  based  on  selected 
lower  degree  functions  and  other  blocks  associated  with  vertices,  edges,  and  faces  of  the 
elements. 

A  bound  on  the  condition  number  independent  of  the  number  of  elements  is  proved, 
and  reasonably  low  growth  of  the  condition  number  has  been  observed  in  practice, 
including  thin  elements,  i.e.,  plates  and  shells.  The  selection  of  the  specific 
preconditioner  adapts  itself  to  the  local  characteristics  of  the  problem  with  no  user 
intervention  required. 


I  MAREK 

Approximations  of  the  principal  eigcnelemenls  of  a  class  of  transport  operators 

Under  some  reasonable  restrictions  on  the  scattering  indicatrix  and  the  cross  sections  a 
class  of  transport  equations  with  boundary  conditions  for  incoming  particles  is 
investigated,  llie  simplest  approximate  spaces  are  proposed  in  which  the  approximating 
operators  preserve  certain  positivity  properties  of  the  operators  approximated. 
Appropriate  error  bounds  are  established. 


A  K  MOHAMMED,  M  H  BALUCH  and  S  T  GOMAA 

Finite  element  modelling  of  deep  beams  using  a  new  refined  theory 

A  recent  refined  theory  that  incorporates  effects  of  transverse  shear,  transverse  normal 
stress  and  transverse  normal  strain  is  used  in  the  finite  element  formulation  for  flexure 
of  deep  beams.  The  finite  element  equations  for  both  beam  bending  and  inplane 
problems  were  developed  using  the  Galerkin  finite  element  method.  Several  problems 
of  uniformly  loaded  deep  beams  with  different  support  conditions  were  solved  to 
evaluate  the  performance  of  the  proposed  finite  element  model.  The  results  are 
compared  with  closed  form  and  elasticity  solutions  whenever  they  exist. 


A  K  MOHAMMED,  M  H  BALUCH  and  S  T  GOMAA 

Finite  element  modelling  of  thick  isotropic  plates  using  a  new  refined  theory. 
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A  recent  refined  theory  that  incorporates  effects  of  transverse  shear,  transverse  normal 
stress  and  transverse  normal  strain  is  used  in  the  finite  element  formulation  of  flexure 
of  isotropic  plates.  The  finite  element  equations  for  both  plate  bending  and  inplane 
problems  were  developed  using  Galerkin  finite  element  methods.  Several  problems  of 
uniformly  loaded  thick  plates  with  different  support  conditions  were  solved  to  evaluate 
the  performance  of  the  proposed  finite  element  model.  The  results  are  compared  with 
Mindlin/Reissner  and  elasticity  solutions  whenever  they  exist 


M  R  MOKHTARZADEH-DEHGHAN,  D  J  STEPHENS  and  L  YU 

Finite  element  analysis  of  turbulent  flow  through  an  orifice  meter  with 
a  concentric  obstruction  of  conical  shape  at  its  centre 

Orifice  meters  provide  a  simple  and  cheap  method  of  measuring  flow  rates.  The 
problem  of  flow  through  an  annular  constriction,  namely  a  simple  orifice  plate  modified 
by  a  concentric  plug,  arises  in  flow  meters  which  are  designed  to  adjust  the  flow  area 
with  variations  in  the  flow  rate.  Such  adjustment  is  done  by  allowing  the  plug  to  move 
in  the  flow  direction  and  therefore  change  the  available  area  through  which  the  fluid 
flows.  In  a  simple  orifice  plate  the  relationship  between  the  flow  rate  and  the  pressure 
difference  across  the  orifice  is  non-linear  and  therefore  the  measured  pressure 
difference,  and  the  accuracy  of  the  meter,  decreases  rapidly  as  the  flow  rate  decreases. 
If  the  flow  area  is  allowed  to  change,  by  careful  design  of  the  plug  profile,  it  is  possible 
to  obtain  a  linear  relationship  between  the  flow  rate  and  the  pressure  drop  and  maintain 
the  accuracy  over  a  wider  range  of  flow  rates.  The  disadvantage  in  this  case  is  the 
increase  in  cost  and  higher  permanent  pressure  loss. 

The  paper  describes  a  finite  element  study  of  turbulent  flow  through  a  simplified 
model,  consisting  of  an  orifice  plate  with  a  concentric  conical  plug  at  its  centre.  The 
mathematical  model  is  based  on  the  time-averaged  equations  of  the  continuity  and 
momentum  in  a  cylindrical  polar  coordinate  system  and  the  k  -  e  turbulence  model.  The 
computer  program  FIDAP  has  been  used  as  the  computational  tool. 

The  velocity  and  pressure  fields,  the  pressure  variation  along  the  pipe  wall  and  the 
pressure  losses  across  the  device  and  the  corresponding  discharge  coefficient  are 
presented  and  discussed.  Comparisons  are  made  with  the  results  obtained  for  a  simple 
orifice  meter. 


P  B  MONK  and  K  PARROTT 

A  dispersion  analysis  of  finite  element  methods  for  Maxwell’s  equations 

A  variety  of  different  finite  element  methods  are  currently  in  use  for  the  solution  of 
Maxwell's  equations.  R  L  Lee  and  N  K  Madsen  (1990)  amongst  others  have  made 
comparisons  between  time  domain  methods  for  these  equations.  An  important  indicator 
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of  the  merit  of  a  given  time  domain  method  is  the  accuracy  with  which  it  propagates 
plane  electromagnetic  waves.  This  is  best  established  through  the  calculation  of  the 
method's  discrete  dispersion  relation.  We  have  carried  this  calculation  out  in  two 
dimensions  on  uniform  triangulations,  consisting  of  either  equilateral  triangles  or 
uniform  right  triangles;  (see  Monk,  1992,  for  similar,  but  less  detailed,  results  on 
quadrilateral  elements). 

On  the  basis  of  this  dispersion  analysis,  it  is  clear  that  some  combinations  of  finite 
element  spaces  are  not  very  attractive  (linear-constant,  edge(l)-Edge(l),  edge(2)- 
constant)  and  some  methods  need  more  analysis  before  becoming  truly  reliable  (non¬ 
conforming-constant).  Two  methods,  the  linear-linear  scheme  and  the  edge(l)-constant 
scheme  have  significant  advantages  and  disadvantages.  The  linear-linear  scheme  has 
excellent  dispersion  properties  on  both  meshes  studied  here  and  can  be  mass  lumped  (at 
the  cost  of  a  decrease  in  accuracy).  However,  it  is  complicated  to  deal  with  boundary 
conditions,  and  the  divergence  conditions  on  the  field  are  only  approximated.  On  the 
other  hand,  the  edge(l)-constant  method  is  more  sensitive  to  the  mesh  showing  fourth 
order  convergence  on  the  equilateral  mesh,  but  only  second  order  on  the  right-angled 
mesb.  Furthermore  the  method  has  unexpected  parasitic  modes.  Mass  Lumping  of  the 
edge-constant  scheme  is  possible  provided  the  mesh  contains  only  strictly  acute  triangles. 
On  the  other  hand  the  method  deals  with  boundary  conditions  and  divergence  conditions 
in  a  natural  way.  These  results  accord  with  the  suggestion  by  G  Mur  (1985,  1992)  that 
optimal  approximations  can  be  obtained  by  selecting  between  these  two  families  of 
elements,  using  appropriate  criteria  reflecting  geometry  and  material  property 
considerations. 

Computation  of  the  dispersion  relations  were  carried  out  using  Mathematics.  These 
results  are  applicable  in  three  dimensions  if  prismatic  elements  are  used  to  discretize 
Maxwell’s 'equations.  The  more  interesting  case  of  tetrahedral  elements  remains  to  be 
analyzed. 
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R  M0CKE  and  J  R  WHITEMAN 

Remarks  on  a-posteriori  error  estimation  for  finite  elasticity 

Several  methods  of  a-posteriori  estimation  of  the  discretization  error  in  finite  element 
solutions  are  well  established  and  widely  used  in  engineering  practice  for  linear  boundary 
value  problems.  In  contrast  we  are  concerned  here  with  finite  elasticity.  The  main 
features  of  finite  elasticity  are  first  summarized.  Using  the  residuals  in  the  equilibrium 
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conditions  (BabuSka  approach)  an  error  measure  in  the  energy  norm  is  proposed  and 
physically  interpreted.  The  effe-'s  arising  from  the  nonlinear  displacement-strain 
mapping  and  the  accuracy  of  the  proposed  error  measure  are  demonstrated  by  a 
numerical  example.  The  estimated  error  distribution  is  then  used  to  control  an  fa- 
adaptive  finite  element  procedure. 


J  NEDOMA 

Finite  element  analysis  in  nuclear  safety 

In  order  to  produce  criteria  for  secure  siting,  location  and  protection  of  operation  of 
high  level  radioactive  waste  repositories  (HLRWRS)  and  nuclear  power  plants  (NPPS) 
mathematical  simulations  must  be  undertaken  of  geodynamic  and  geomechanic 
processess  in  both  the  global  and  the  local  areas,  where  HLRWRS  or  NPPS,  RESP,  will 
be  situated.  In  the  contribution  the  finite  element  analyses  of  mathematical  models 
based  on  the  contact  problems  of  Signorini  type  in  thermoelasticity  and  on  the 
nonstationary  incompressible  thermo-Bingham  problem  will  be  discussed.  The  methods 
for  monitoring  the  stress-strain  field  in  situ,  of  the  determination  of  the  optimal  design, 
of  the  measuring  tensometric  probes  and  of  the  postprocessing  of  measured  data  based 
on  the  FEM  will  also  be  briefly  discussed.  Existence  and  convergence  theorems  will  be 
proved,  and  the  algorithms  will  be  presented. 


J  NEDOMA 

FEM  analysis  of  artificial  substitutes  of  human  joints  and  their  optiml  design 

The  actual  situation  in  human  joints  under  loading  is  much  more  complicated  than  that 
described  with  classical  methods  of  biomechanics  such  as  Bombclli  (1983).  Some  new 
ideas  of  biomechanics  of  human  joints  and  of  their  artificial  substitutes  were  given  in 
Nedoma  and  Stehlik  (1989),  Nedoma  (1991).  These  are  based  on  contact  problems  in 
elasticity  and  thermoelasticity  (plasticity). 

In  the  contribution  mechanical  processes  which  take  place  during  static  burdening  in 
the  contact  area  of  the  acetabulum  and  the  head  of  the  hip  and  of  their  artificial 
substitutes  -  the  artificial  acetabulum  and  the  head  of  the  endoprosthesis  -  will  be 
discussed.  The  mathematical  problem  introduced  in  the  contribution  represents  a 
simulation  of  the  function  of  human  joints  and  their  total  substitutes  by  total 
endoprostheses  (TEP).  The  model  problems  lead  to  coercive  and  semicoe reive  contact 
problems  in  quasi-coupled  thermoelasticity.  The  FEM  is  used  for  numerical  analyses  of 
such  model  problems.  Hie  optimal  design  of  the  TEPs  based  on  the  FEM  analyses  of 
the  contact  problems  in  thermoelasticity  is  also  briefly  discussed  and  existence  and 
uniqueness  of  the  solutions  of  the  model  problems  and  a  convergence  theorem  will  be 
proved. 
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P  NEITTAANMAKI  and  V  RIVKIND 

Finite  element  method  for  model  problems  of  changing  phase 

We  consider  the  model  of  multiphase  flow  described  through  full  coupled  Navier-Stokes 
and  Stephan  equations  with  special  boundary  conditions.  Proof  of  the  existence  and  the 
uniqueness  is  given  (with  an  assumption  on  the  limitation  of  initial  parameters  of  the 
problem).  A  finite  element  method  for  the  numerical  solution  is  given.  The  existence 
and  uniqueness  of  the  approximate  solutions  is  proved  (with  assumptions  of  the 
uniqueness  and  smoothness  of  the  exact  solution).  Proofs  of  convergence  and  of  the  rate 
of  convergence  of  the  approximate  solution  to  the  exact  one  are  given.  Several  iterative 
methods  such  as  the  full  and  partial  linearization,  Newton’s  method,  and  the  domain 
decomposition  method  are  proposed.  The  case  of  many  media,  including  media  with 
large  relations  between  parameters  of  sub-media,  is  discussed. 

I 

H  ORS 

Numerical  simulation  of  flow  circulation  in  the  golden  horn 

A  numerical  model  of  the  flow  in  the  Golden  Horn  is  presented.  The  model  is  based 
on  shallow  water  theory.  Considering  the  complex  geometry  of  the  domain  of  interest, 
the  governing  set  of  partial  differential  equations,  with  appropriate  boundary  conditions, 
is  solved  by  the  finite  element  method.  Given  the  inner  sea  nature  of  the  Marmara  Sea 
with  minimal  tidal  effects,  the  steady  state  case  is  considered.  The  computational 
domain  is  discretized  into  140  triangular  elements  with  111  node  points.  Different 
interpolation  functions  are  tried  in  the  context  of  the  Galerkin  approach.  The  final 
system  of  equations  is  solved  by  a  Cholesky  factorization  method.  Flow  field  values  thus 
obtained  are  compared  with  on  site  measurement  values.  The  ultimate  goal  of  the 
project,  supported  by  Bogazifi  University  Research  Fund,  being  to  study  the  pollution 
mechanisms  of  the  Golden  Horn  and  ways  of  cleaning  it,  a  second  code  is  written  to 
compute  the  dispersion  rates  of  the  pollutants.  Values  of  the  velocity  field  previously 
obtained  are  used  in  incorporating  convection  terms.  Different  wind  force  distributions 
and  pollution  sources  are  also  tried.  Computations  are  carried  out  on  a  HP  720 
workstation  and  results  are  presented  graphically. 
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F  PENZEL 

Error  estimates  for  discretized  boundary  element  methods  for 
three-dimensional  crack  problems 

The  displacement  in  an  unbounded  elastic  solid  with  a  given  internal  crack  under  normal 
traction  can  be  calculated  with  the  aid  of  boundary  element  methods. 

Let  us  assume  that  the  crack  surface  T  C  {(xpcjO)  :  (XpX2)  €  M*}  is  a  polygon.  The 
crack  opening  displacement  u  £  [ht'(T)}>  is  known  to  be  a  solution  of  the  hyper¬ 
singular  boundary  integral  equation  Du  =  v,  v  6  is  the  known  traction  vector 

and  D  is  a  well-known  pseudodifferential  operator  of  order  + 1. 

The  COD  u  is  approximated  by  a  Galerkin  scheme  using  a  regular  triangulation  of  T  and 
piecewise-linear  finite  elements.  Let  h  be  the  supremum  of  the  diameters  of  the 
elements  in  a  fixed  triangulation.  If  the  traction  data  v  is  a  sufficiently  smooth  function, 
then  it  is  known  that  the  solution  uk  of  the  Galerkin  scheme  converges  to  u  with  an 
order  of  A**’*,  for  j/  >  0  chosen  arbitrarily  small. 

Here  we  present  a  fully  discretized  Galerkin  scheme  obtained  by  use  of  a  quadrature 
formula  in  the  definition  of  the  scalar  products.  We  reduce  the  analysis  of  convergence 
of  the  fully  discretized  Galerkin  method  for  the  hypcrsingular  integral  equation  Du  * 
v  to  the  analysis  of  the  weakly  singular  integral  operator  -1.  Under  the  assumption  that 
the  quadrature  error  is  sufficiently  small  we  prove  convergence  of  order  of  the 
solution  of  the  fully  discretized  Galerkin  method.  By  use  of  explicitly  determined  error 
bounds  we  formulate  restrictions  on  the  quadrature  formulas  which  are  sufficient  to 
ensure  this  convergence  result. 


E  PRIOLO  and  G  SERIANI 

Low-  and  high-order  FEM  in  seismic  modelling:  Experience  and  perspectives. 

In  this  paper  we  describe  some  experience  in  improving  the  computational  efficiency  of 
a  finite  element  code  based  on  a  global  approach,  used  for  seismic  modelling  in  the 
framework  of  geophysical  oil  exploration.  Applications  in  this  area  are  characterized  by 
models  that  have  complex  geometries  and  heterogeneous  structures,  and  can  be 
successfully  handled  by  the  finite  element  method.  Boundary  conditions  such  as  a  free 
surface  can  easily  be  taken  into  account.  The  main  drawbacks  of  a  classical  approach 
based  on  the  assembly  of  global  matrices  are  the  low  computational  efficiency  and  the 
possible  appearance  cl  spurious  effects. 

Results  coming  from  runs  of  different  methods  and  models  on  a  mini  super 
workstation  APOLLO  DN 10000  are  reported  and  compared.  Low-  and  high-order 
elements  are  used  to  solve  both  the  2D-acoustic  and  2D-elastic  wave  equations.  Tune 
integration  is  performed  with  implicit  Newmark  two-step  method.  The  use  of 
compressed  storage  is  shown  to  greatly  reduce  not  only  memory  requirement  but  also 
computing  time.  With  Chebyshev  spectral  elements  great  accuracy  can  be  reached  with 
almost  no  numerical  artifacts.  For  spectral  elements  the  static  condensation  of  internal 
nodes  reduces  memory  requirements  and  CPU  time,  and  is  well  suited  to  parallel 
computing  in  a  distributed  environment.  With  this  technique  the  solution  of  the  linear 
system  is  performed  with  a  hybrid  method  by  using  a  preconditioned  conjugate  gradient 
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for  tbe  Sbur  complement  and  a  direct  method  for  the  internal  node  matrices.  A  speed¬ 
up  of  about  3  and  a  saving  in  memory  of  about  25%  is  reached  in  this  way.  The  matrix- 
by-vector  product  of  both  sparse  and  full  matrices  is  a  fundamental  operation  for  which 
an  efficient  programming  and  storage  technique  is  decisive. 

The  size  of  2D-geological  models  that  can  be  handled  in  a  reasonable  time  on 
computers  of  this  kind  is  nowadays  hardly  sufficient,  whereas  3D  modelling  is  completely 
unfeasible.  However,  the  introduction  of  new  FEM  algorithms  (Legendre  spectral 
elements,  element-by-element  approach,  and  domain  decomposition)  coupled  to  the  use 
of  new  emerging  computer  architectures  are  encouraging  for  the  future. 


E  D  PROVIDAS 

A  simple  and  efficient  solution  for  the  facet  approximation  of  shells 

A  method  for  adding  drilling  rotation  degrees  of  freedom  to  the  constant  strain  triangle 
is  presented.  The  method  is  free  from  any  deficiencies  and  does  not  alter  the  strain 
energy  of  the  constant  strain  triangle.  Superposition  with  the  displacement  version  of 
the  hybrid  bending  element  yields  tbe  required  six  degrees  of  freedom  per  vertex. 
Several  standard  shell  problems  are  solved.  Attention  is  focussed  on  testing  the  capacity 
of  the  element  to  simulate  inextensional  bending  of  real  arbitrary  shells.  For  this 
purpose  a  displacement  prescribed  patch  test  and  a  matrix  procedure  are  employed.  The 
element,  in  contrast  to  other  similar  elements,  is  proven  to  respond  satisfactorily  to  this 
type  of  deformation  modes. 

i 

W  RACHOWICZ 

Anisotropic  h-adaptive  finite  element  method  for  hypersonic  viscous  flows 
with  shock-boundary  layer  interaction 

Viscous  hypersonic  flows  are  analysed  using  second  order  Taylor-Galerkin  scheme.  An 
h-adaptive  version  of  the  finite  element  method  is  used  to  approximate  the  solutions. 
Tbe  method  supports  anisotropic  b-refinements,  i.e,  tbe  possibility  of  breaking  of 
quadrilateral  elements  in  one  of  the  two  possible  directions.  The  strategy  of  selecting 
elements  to  be  refined  and  the  direction  of  the  refinement  is  based  on  minimization  of 
the  interpolation  error  and  it  involves  estimation  of  second  order  derivatives  of  the 
solution.  The  anisotropic  strategy  is  especially  effective  in  boundary  layers  where 
solutions  are  almost  one-dimensional.  The  method  allows  for  accurate  resolution  of 
aerothermal  loads  with  considerably  smaller  computational  effort  than  with  the  standard 
adaptivity.  Numerical  simulation  of  a  shock-boundary  layer  interaction  problem 
illustrates  the  approach. 


55 


A  RAMAGE  aid  A  I  WATHEN 

On  preconditioning  for  finite  element  equations  on  irregular  grids 

One  efficient  algorithm  for  solving  Galerkin  finite  element  equations  is  the 
preconditioned  conjugate  gradient  method.  A  large  number  of  preconditioning  strategies 
have  been  proposed  and  these  are  generally  analysed  and  compared  using  model 
problems:  simple  discretisations  of  Laplacian  operators  on  regular  computational  grids, 
generally  in  two  space  dimensions.  Many  of  these  techniques  have  attractive  theoretical 
convergence  estimates  on  such  model  problems  and  it  is  principally  for  geometrically 
irregular  (non-model)  problems  that  the  applicability  and  economy  of  preconditioned 
conjugate  gradient  methods  is  most  useful.  This  is  particularly  true  for  problems  on 
irregular  unstructured  three-dimensional  grids. 

Developing  rigorous  theory  for  irregular  problems  has  proved  to  be  more  difficult  than 
in  analogous  regular  cases  and,  as  a  result,  many  aspects  of  working  with  unstuctured 
finite  element  grids  are  a  lot  less  well  understood.  Here  we  present  some  theoretical 
results  which  apply  to  any  finite  element  grid  regardless  of  irregularity.  We  extend  the 
work  of  the  second  author  to  find  easily  computed  eigenvalue  bounds  for  a  class  of  finite 
element  matrices  on  irregular  grids.  In  addition,  we  obtain  (weak)  bounds  on  the  interior 
eigenvalues  of  a  general  finite  element  matrix.  These  are  useful  in  the  analysis  of  various 
preconditioners  for  irregular  finite  element  problems.  We  also  present  a  result 
concerning  the  effect  of  a  certain  class  of  element-based  preconditioners  on  the  matrix 
condition  number,  which  is  again  relevant  to  the  rate  of  PCG  convergence. 


A  RANJBARAN 

A  precis  of  developments  in  modelling  embedded  reinforcement  for  the 
finite  element  analysis  of  reinforced  concrete  structures. 

To  date  three  alternatives  are  available  for  modelling  of  reinforcement  in  reinforced 
concrete  elements,  i.e.  discrete,  smeared  and  embedded. 

In  this  paper  various  techniques  of  modelling  reinforcement  are  critically  reviewed  and 
their  merits  and  shortcomings  are  highlighted.  The  embedded  models  are  put  into  two 
general  categories,  i.e.  locally  embedded  model  (LEM)  and  globally  embedded  model 
(GEM).  When  the  trajectory  of  a  reinforcement  in  isoparametric  coordinates  of 
concrete  is  assumed  as  a  priori,  the  model  is  called  LEM.  On  the  other  hand,  if  the 
trajectory  in  global  coordinates  is  known  the  term  GEM  is  used.  A  totally  general  GEM 
is  proposed.  The  relevant  mathematical  formulation  for  computation  of  stiffness  and 
load  contribution  of  a  reinforcing  bar  to  those  of  concrete  are  derived.  The  two-  and 
three-dimensional  problems  are  treated  in  a  unified  manner.  Moreover,  it  is  shown  that 
all  previous  embedded  models,  available  in  the  literature,  can  be  considered  as  special 
cases  of  the  proposed  model.  The  proposed  model  is  implemented  in  a  finite  element 
program.  Numerical  examples  are  included  to  verify  the  efficacy  of  the  model. 
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E  RANK 

A  zooming-technique  using  a  hierarchical  hp-version  of  the  finite  element  method 

The  hp-version  of  the  finite  element  method  combines  local  mesh  refinement  with  an 
increase  of  the  polynomial  order  of  the  shape  functions.  It  has  been  shown  theoretically 
and  in  many  numerical  examples,  that  exponential  rate  of  convergence  in  energy  norm 
can  be  obtained  for  linear  elliptic  boundary  value  problems.  The  bp-version  has  also 
been  applied  successfully  to  more  general  problems  like  reaction  diffusion  equations  or 
nonlinear  Navier-Stokes  equations.  Recently  a  variant  of  the  hp-version  as  a 
combination  of  a  high  order  approximation  with  a  domain  decomposition  method  has 
been  suggested  [1]  by  the  author.  This  ’hp-version’  is  very  similar  to  Fish’s  ’s-version’ 
[2|  using  a  superposition  of  independent  finite  element  meshes. 

The  basic  idea  can  be  explained  as  follows.  In  a  first  step  of  the  analysis  a  pure  p- 
version  approximation  is  performed  on  a  coarse  finite  element  mesh.  Controlled  by  user 
interaction  or  by  an  a  posteriori  error  estimation  the  coarse  mesh  is  then  covered 
partially  by  a  geometrically  independent  fine  mesh.  On  this  second  mcsb  a  low  order 
approximation  is  performed  and  the  global  approximation  is  defined  as  the  hierarchical 
sum  of  the  p-approximation  on  the  coarse  mesh  and  the  h-approximation  on  the  fine 
mesh.  Global  continuity  of  the  finite  element  solution  can  be  guaranteed  by  imposing 
homogeneous  conditions  at  the  fine  mesh  boundary.  The  hierarchical  nature  of  the 
approximation  also  reflects  in  the  structure  of  the  arising  linear  equation  system  and  can 
be  used  in  an  efficient  solution  algorithm.  As  an  example  a  composed  mesh  and  contour 
lines  for  a  potential  flow  around  a  well  will  be  presented.  On  the  coarse  mesh 
polynomial  degree  p  =  4  is  used,  the  fine  mesh  approximation  is  performed  with  linear 
elements. ' 

The  paper  discusses  algorithmic  details  of  the  suggested  approach  and  addresses  the 
question  of  how  this  domain  decomposition  can  be  implemented  without  major 
modifications  into  existing  finite  element  codes.  An  extension  of  the  method  offers  also 
the  possibility  of  'zooming'  areas  of  critical  solution  behaviour.  It  is  possible  to  compute 
the  global  p-version  approximation  with  'averaged'  material  parameters,  and  only  locally 
to  resolve  the  structure  in  detail.  It  will  be  shown  how  the  hierarchical  approach  offers 
the  possibility  of  a  consistent  modelling  of  such  local-global  behaviour.  In  several 
numerical  examples  the  efficiency  and  accuracy  of  the  method  will  be  demonstrated. 
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H-G  ROOS 

Nonconforming  exponentially  fitted  elements 

We  consider  the  singularly  perturbed  boundary  value  problem 
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-  cAk  +  bVu  +  cm  =/  in  Q  m 

u  -  0  on  Q 

under  the  assumptions  c  •  112  divb  >  0,  |  b  |  *  0.  There  are  a  lot  of  special  finite 
element  techniques  to  discretize  (1)  in  the  convection  dominated  case:  upwinding 
(Griffiths,  Mitchell),  symmetrization  (Morton),  streamline  diffusion  (Hughes).  But  only 
based  on  some  exponential  fitting  is  it  possible  to  show  uniform  convergence  with  respect 
to  the  singular  perturbation  parameter.  A  first  result  in  this  direction  is  due  to 
O’Riordan  and  Stynes  using  tensor  products  of  usual  lD-L-splines  for  both,  the  trial  and 
the  test  space.  Under  the  restrictive  assumption  b  =  (b,b3)  and  b  =  (b,(x\b3(y)  on 
a  rectangular  domain  they  proved  uniform  0(hja)  convergence. 

For  the  general  case  b  =  (b,(x,y),b/x,y))  it  seems  to  be  adequate  to  use  exponentially 
fitted  splines  which  are  no  longer  globally  continuous.  Therefore  it  is  necessary  to 
investigate  the  corresponding  nonconforming  Petrov-Galerkin  methods. 

First  we  will  present  some  abstract  convergence  theorems  for  nonconforming  Petrov- 
Galerkin  methods  which  generalize  the  well  known  theorems  of  Strang. 

Further  we  will  show  that  these  abstract  results  are  a  suitable  tool  to  handle 
nonconforming  exponentially  fitted  elements.  We  study  the  model  problem 

-at"  +  b(x)u  '  *  c(x)u  m  f(x) 
u(0)  -  u(l)  -  0 

and  choose  the  following  combinations  for  the  trial  and  the  test  spaces: 

(i)  conforming  linear  elements/approximated  (nonconforming)  L ’-splines 

(ii)  approximated  L-splines/approximated  L ’-splines. 

In  both  cases  we  are  able  to  prove  that  the  discrete  problems  admit  unique  solutions. 
The  key  of  the  error  analysis  consists  in  establishing  a  stability  condition 

3  V*  e  Xt  :  *>*(«*,¥>)  ^  mt|uj  |y>|  (mk  >  0) 
for  the  corresponding  bilinear  form 

vmo  -  g1  /;*■  «»♦'  -  fa*  - *  (c,  -  ')■*♦*. 

Finally  we  prove  the  optimal  error  estimations 

(i)  1|  u'  -  u„  |  S  Chw 

(ii)  |  u  ■  uk  H  s  CTi* 

(C  is  independent  of  e)  with  respect  to  the  adapted  energy  norm. 


S  RIPPA  and  B  SCH1FF 


Improving  ihe  accuracy  of  finite-element  solutions  to  two-dimensional  elliptic  problems 

We  present  a  method  for  improving  the  accuracy  of  a  finite-element  solution  of  a  two- 
dimensional  elliptic  problem  over  a  triangular  mesh  by  modifying  the  triangulation  whilst 
keeping  the  mesh  points  unaltered.  The  triangulation  is  modified  by  a  process  of  "local 
optimisation"  in  which  each  quadrilateral  formed  by  two  triangles  sharing  a  common 
edge  is  considered  in  turn.  The  value  of  the  energy  integral  over  the  two  triangles  is 
.compared  with  the  integral  over  the  two  triangles  which  would  be  obtained  had  the  given 
quadrilateral  been  divided  by  the  opposite  diagonal  If  the  latter  integral  is  less  than  the 
former,  we  modify  the  triangulation  by  replacing  the  diagonal  of  the  given  quadrilateral 
by  the  oposite  diagonal.  As  the  remainder  of  the  triangulation  is  unaltered,  the  value 
of  the  energy  integral  over  the  whole  domain  has  decreased.  After  completing  a  sweep 
of  local  optimisation  over  all  of  the  quadrilaterals,  the  finite  element  equations  may  be 
solved  again  over  the  new  triangulation  to  obtain  new  nodal  values.  Alternatively,  before 
re-solving,  additional  sweeps  of  local  optimisation  can  be  carried  out  using  the  same 
nodal  values.  In  either  case,  once  the  equations  have  been  solved  over  the  new 
triangulation,  local  optimisation  is  found  to  lead  to  no  significant  improvement  in  the 
solution.  A  third  possibility  is  to  terminate  the  procedure  without  re-solving  at  all. 

The  method  has  been  tested  for  the  Poisson  equation  using  linear  elements  over  each 
triangle,  and  for  the  biharmonic  equation  using  (C1)  cubic  macroelements  of  Clough- 
Tocher  and  reduced  Clough-Tocher  (normal  derivative  linear  over  the  external  edges  of 
the  macro-element)  type.  Various  problems  with  known  solutions  were  solved,  and  the 
process  often  led  to  considerable  improvement  in  the  results,  particularly  when  the 
solution  has  a  marked  directional  dependence.  For  the  biharmonic  equation,  two 
problems  with  a  boundary  singularity  were  also  solved,  the  determination  of  the  Airy 
stress-function  for  a  rectangular  elastic  plate  containing  an  edge  crack  under  uniform 
tension,  and  the  Stokes  flow  in  a  square  domain  with  one  of  the  four  walls  moving  with 
constant  velocity  (the  "driven  cavity"  problem).  Comparison  with  existing  solutions 
obtained  for  these  two  problems  by  a  variety  of  other  methods  shows  that  the  "local 
optimisation"  technique  results  again  in  a  considerable  improvement  of  the  original 
finite-element  solution. 


C  SCHWAB  and  M  SURI 

The  numerical  resolution  of  boundary  layers  and  locking  effects 
in  the  Reissner-Mindlin  plate  model 

The  Reissner-Mindlin  plate  model  is  characterized  by  both  boundary  layers  and  locking 
as  the  thickness  t  tends  to  zero.  The  boundary  layers  have  the  form  exp{-Xp/t),  where 
A  >  0  is  a  known  constant,  and  p  is  the  distance  to  the  boundary.  For  t  small,  the 
approximation  of  such  Ikayers  can  be  quite  unsatisfactory  unless  the  discretization 
parameter  is  sufficiently  small.  Locking  effects  can  also  cause  the  approximations  to 
deteriorate  when  t  is  small.  These  occur  due  to  the  Kirchhoff  constraint  being  imposed 
on  both  the  true  and  the  approximate  solutions  in  the  limiting  case,  as  t  tends  to  zero. 
In  order  to  accurately  approximate  such  models  numerically,  it  is  essential  that  the 
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numerical  method  under  consideration  be  designed  in  such  a  way  that  both  these  effects 
are  resolved  uniformlym,  i.e.  with  errors  independent  of  the  thickness. 

In  this  talk,  we  present  various  results  for  the  approximation  of  boundary  layers  by  the 
h,  p  and  h-p  finite  element  method.  For  the  h  version  we  show  how  exponentially 
graded  meshes  towards  the  boundary  yield  an  optimal  convergence  rate  independent  of 
t.  We  also  show  that  by  a  suitable  combination  of  mesh  refinement  and  degree  selection 
(the  h-p  method),  one  can  obtain  uniform  exponential  rates  of  convergence. 

Next,  we  discuss  the  locking  and  robustness  properties  of  various  h  and  p  type  finite 
element  methods.  We  quantify  the  exact  degree  of  locking  that  can  be  expected  if 
various  t  riangular  and  rectangular  meshes  are  used.  We  show  that  the p  version  is  free 
of  locking  for  this  problem. 


K  SFGETH 

Numerical  experience  with  grid  adjustment  based  on  a  posteriori  error  estimators. 

Recently,  a  variety  of  techniques  for  a  posteriori  error  estimation  have  been  theoretically 
developed  and  practically  applied.  A  posteriori  estimates  can  serve  as  a  means  for  grid 
adjustment  ensuring  the  optimal  number  and  optimal  distribution  of  grid  points  in  the 
finite  element  method.  We,  however,  need  a  solution  of  the  problem  on  some  grid  to 
construct  a  new,  optimal  grid.  This  approach,  therefore,  is  very  suitable  e.g.  for  solving 
parabolic  partial  differential  equations  by  the  method  of  lines.  The  analysis  of  the 
approximate  solution  at  a  fixed  time  level  then  yields  a  new  grid  to  be  used  for  the  time 
step  leading  to  the  next  time  level. 

In  the  paper,  some  approaches  using  the  idea  of  monitor  equidistribution,  with 
monitors  based  on  a  posteriori  error  bounds  as  well  as  with  those  independent  of  error 
estimators,  are  presented  together  with  numerical  experience. 


D  SILVESTER  and  A  J  WATHEN 

Efficient  preconditioners  for  incompressible  flow 

Having  an  efficient  Stokes  solver  is  an  important  requisite  of  many  Navier-Stokes 
solution  algorithms.  The  Stokes  operator  is  self-adjoint  so  iterative  solution  methods  are 
an  attractive  approach  for  large  problems.  The  only  impediment  to  efficiency  is  the 
indefinite  nature  of  the  discrete  system. 

In  this  talk,  some  new  theoretical  results  will  be  presented  showing  the  effect  of 
preconditioning  discrete  Stokes  systems.  Our  analysis  covers  a  variety  of  preconditioners, 
including  simple  diagonal  scaling,  incomplete  factorisations  and  multilevel 
preconditioning.  Our  results  apply  to  stable  approximations  of  the  Stokes  operator, 
typically  based  on  staggered  grids,  and  also  to  recently  developed  stabilised 
approximations:  equal  order  finite  element  methods  or  finite  difference  methods  on 
unstaggered  grids. 
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M  SLODlCKA 


Numerical  solution  of  a  parabolic  equation  with  a  weakly  singular 
positive  type  memory  term 

There  exist  many  models  from  heat  conduction  in  materials  with  memory,  population 
dynamics  and  visco-elasticity,  which  can  be  described  by  integro-differential  equations. 
Usually,  the  memory  term  depends  on  the  solution  u  or  on  its  spatial  derivatives.  But 
there  exist  cases  where  this  memory  term  can  depend  on  the  time  derivative  of  u. 
Integro-differential  equations  of  this  nature  appear  in  describing  diffusion  models  for 
fractured  media  (c.f.  Hornung-Showalter  [1]). 

We  discuss  the  numerical  solution  of  the  following  problem 

ii,(0  -  An(f)  +  f'0(t-s)up)d3  0)  to  Q.  t  >  0  , 

n  -  0  ob  30  ,  t  >  0  , 
n(0)  -  v  in  O  , 


where  {lisa  bounded  convex  polyhedral  domain  inW'fd  Z  1),  v  e  fl'fD).  Our  integral 
kernel  (in  the  memory  term  is  supposed  to  be  weaking  singular  and  positive  type. 

We  use  the  backward  Euler  method  for  discretization  in  time  and  the  Galerkin  finite 
element  method  for  discretization  in  space.  We  prove  the  convergence  of  our 
approximation  scheme  in  some  functional  spaces,  the  existence  and  uniqueness  of  the 
solution,  and  produce  some  regularity  results  for  exact  solution. 
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E  STEIN  and  F-J  BARTHOLD 

Design  optimization  of  rubber  elastic  structures 

Sensitivity  analysis  of  physical  and  geometrical  nonlinear  behaviour  is  an  outstanding 
research  field  within  design  optimization.  In  this  lecture  sensitivity  analysis  of  rubber 
elastic  structures  and  its  efficient  implementaiton  into  a  design  optimization  procedure 
is  outlined. 

Based  on  continuum  mechanics  at  finite  strains,  essential  remarks  on  the  finite 
element  method  for  nearly  incompressible  rubber  elastic  materials  and  on  the  efficient 
use  of  adaptive  mesh  refinement  are  discussed.  Then  the  sensitivity  analysis  for  rubber 
elastic  materials  with  respect  to  changes  of  geometry  and  material  parameters  is 
described  using  a  parametric  model  consisting  of  space  coordinates,  time  and  design 
variables.  Due  to  large  strains  occurring  in  technical  applications  of  rubber  materials  a 
current  configuration  approach  is  used  both  in  structural  and  sensitivity  analysis. 

The  concept  of  a  design  optimization  procedure  and  integrated  algorithms  for 
Computer  Aided  Design  (CAD),  Finite  Element  Method  (FEM),  Mathematical 
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Programming  and  Sensitivity  Analysis  is  explained.  A  highly  modular  subprogram  and 
data  base  concept  for  these  methods  is  used  to  implement  efficiently  the  analytical 
derived  sensitivities  into  the  FE  research  program  INA-OPT  (INelastic  Analysis  and 
OPTimization)  developed  at  our  institute.  Applications  of  the  optimization  procedure 
to  material  parameter  identification  problems  in  tyre  layout  are  shown. 

The  methodology  for  the  approximate  solution  of  the  inverse  problem  was  used  for 
the  identification  of  material  parameters  in  the  Ogden  material  equation  which  can  also 
be  understood  as  a  validation  of  the  constitutive  model  for  the  tyre  material  investigated. 

An  optimal  design  of  steel  cords  in  a  tyre  layout  based  on  a  minimal  internal  energy 
equilibrium  state  is  derived  using  the  optimization  procedure  described. 


R  STENBERG 

On  linear  and  bilinear  elements  for  Reissner-Mindlin  plates 
We  consider  the  following  elements: 

1.  The  MITC4  analyzed  by  Brezzi  and  Bathe  in  the  MAFELAP  1984  Conference  for 
the  case  of  rectangular  elements.  We  give  the  error  analysis  for  the  more  general  class 
of  meshes  for  which  the  "Q,-Pa"  Stokes  element  was  analyzed  by  Pitkhranta  and  Stenberg 
in  the  same  conference. 

2.  We  consider  the  linear  triangular  MITC  element  in  which  the  deflection  and  rotation 
are  linear  and  the  reduction  operator  is  given  by  the  lowest  order  rotated  Raviart- 
Tbomas  element.  We  give  an  error  analysis  for  the  case  when  the  triangular  mesh  is 
obtained  from  a  quadrilateria!  one  by  drawing  the  two  diagonals  in  each  quadrilateral. 
The  error  analysis  can  now  be  performed  for  the  same  quadrilateral  meshes  as  for  the 
M1TC4.  Our  analysis  seems  to  explain  the  good  behaviour  reported  by  Hughes  and 
Taylor  in  MAFELAP  1981  for  a  very  similar  element  As  a  byproduct  we  also  give  an 
error  analysis  of  the  well  known  linear  constant  strain  triangle  with  a  "criss-cross"  mesh 
for  an  incompressible  material. 

3.  We  show  that  both  of  the  above  elements  can  be  modified  so  that  they  are  optimally 
convergent  without  any  restrictions  on  the  meshes. 


A  TESSLER,  H  R  RIGGS  and  S  C  MACY 

Application  of  a  variational  method  for  computing  smooth  stresses,  stress  gradients 
and  error  estimation  in  finite  element  analysis 

The  displacement  finite  element  method  is  known  to  be  deficient  as  far  as  the  accuracy 
of  strain  and  stress  predictions  is  concerned.  This  aspect  is  particularly  regrettable  since 
accurate  predictions  of  strains  and  stresses  are  necessarily  needed  in  structural  design 
and  failure  predictions.  Moreover,  in  the  process  of  adaptively  refiniog  the  finite 
element  discretization  to  achieve  converged  solutions  in  an  automated  and  systematic 
way,  a-posteriori  estimation  of  errors  is  required.  Error  estimation  is  commonly  pursued 
by  way  of  computing  improved  stress  solutions  which  are  used  to  measure  the  goodness 
of  the  finite  element  solution. 

The  source  of  the  computational  difficulty  in  recovering  strains  and  stresses  is  well 
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established.  The  displacement  field  is  approximated  with  piecewise  C°-continuous  shape 
functions  ensuring  continuity  of  displacements  along  element  boundaries.  The  strains 
(and  stresses)  are  the  dependent  quantities  that  are  generally  represented  by  first 
gradients  of  the  displacements;  hence,  they  are  only  C*  continuous,  exhibiting 
nonphysical  discontinuity  along  element  boundaries.  Reasonably  accurate  strain  (stress) 
results  can  generally  be  recovered  at  the  limited  number  of  discrete  locations  within  a 
finite  element,  these  being  optimal  (or  Barlow)  points.  Since  these  points  are  in  the 
element  interior,  and  for  the  lowest-order  elements  there  is  only  a  single  centroidal 
point,  reliable  stgrain  (stress)  values  on  the  element  boundaries  are  generally  not 
available. 

Recently,  Tessler  and  co-workers  developed  a  smoothing  formulation  based  on  a 
variational  principle  which  combines  weighted  least-squares  and  penalty-constraint 
functionals  in  a  single  variational  form  and  uses,  in  the  general  three-dimensional  case, 
four  independent  field  variables  -  the  smoothed  quantity  (strain  or  stress)  and  its  three 
orthogonal  gradients.  The  unique  characteristic  of  this  approach  is  that  the  resulting 
stress  field  is  globally  C'-continuous. 

The  mathematical  formulation  can  be  briefly  stated  as  foltows: 

Le»  d,  -  dfcr,)  (x,  -  (x,‘.x’.x,,)l  p  -  1X-...N)  represent  a  set  of  discrete  data,  such 

as  Barlow-point  strains  (stresses),  and  let  fl  =  (x  =  (x',x\x5)  C  SH’J  denote  the  domain 
of  the  finite  element  discretization  defined  in  the  three-dimensional  orthogonal  frame. 
We  seek  a  smooth,  scalar-valued  continuous  (stress)  function  o(x)  which  can  best 
represent  the  discrete  data  by  casting  the  problem  as  a  minimization  of  the  following 
least-squares/penalty-constraint  functional 

-  0 


with  ♦  given  as 


•  -  E  [VBM*  +  k  /0l(° i-o,)2  ♦  <«,-«>*  *  (,) 

where  N  is  the  total  number  of  data  {mints,  xp  denotes  the  position  vector  of  the  data 
point,  A  is  a  large  penalty  number,  0,,  02,  and  0 3  are  independent  continuous  functions 
also  defined  on  O,  and  o{  =  d/dx,  a  (i  =  1,2,3).  The  minimization  of  4>  is  performed 
with  respect  to  the  coefficients  in  cr  and  B,  which  serve  as  the  unknowns  in  this  problem. 
The  first  term  in  (1)  represents  a  discrete  least-squares  functional)  with  the  term 

«,  *  0,  -  oty  (2) 

representing  the  error  in  op  as  compared  with  the  smooth  solution  o(xr).  The  second 
term  in  (1)  is  a  penalty  constraint  functional  which,  for  A  -*  »,  yields  the  following 
constraint  relations 

•  (I-1A3)  on  0  -  (3) 

This  means  that,  for  all  practical  purposes  as  long  as  A  is  very  large,  0,  (i  =  1,2,3) 
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represent  the  first  order  gradients  of  a  with  respect  to  the  orthogonal  coordinates. 

The  present  paper  revists  the  smoothing  variational  formulation  for  the  purpose  of 
demonstrating  its  effectiveness  and  robustness  in  the  three  applicaiton  areas:  (1) 
recovery  of  accurate,  globally  C1  continuous  stresses  (2)  computations  of  ply-level 
transverse  shear  stresses  in  laminated  composite  shell  analysis,  and  (3)  computations  of 
reliable,  exact  local-equilibrium  error  estimators  needed  in  adaptive  mesh  refinement 
analyses.  The  success  of  the  analyses  (2)  and  (3)  is  related  directly  to  the  effectiveness 
of  the  smoothing  method  to  recover  stress  (strain)  gradients  from  a  conventional, 
displacement-based  finite  element  analysis. 

J  L  TORRES 

Finite  element  solution  of  unsteady-state  mass  transfer  through  a  stationary  liquid 

Mass  transfer  by  diffusion  of  one  liquid  component  through  another  can  be  modelled  by 
the  expression: 

Sc 

-  V .  cDAt'VyA  *  V .  eAv  +  ~  -  JtA  ■  0  .  0) 

If  constant  liquid  density  and  constant  diffusivity  are  assumed,  and  if  the  diffusion  occurs 
through  a  stationary  layer,  the  expression  reduces  to  the  so-called  Fick’s  Second  Law  of 
Diffusion  for  non-reacting  systems: 

is  -  ■  (2) 

The  calculation  of  a  concentration  profile  as  a  function  of  time  and  position  in  such 
systems  is  usually  carried  out  using  equation  (2),  mainly  because  the  rigorous  application 
of  (1)  can  lead  to  computational  complications. 

The  assumptions  that  lead  to  equation  (2)  are  usually  appropriate,  with  one  important 
exception.  When  the  diffusion  takes  place  in  a  stationary  liquid  layer,  in  many  cases  the 
diffusivity  is  a  strong  function  of  concentration,  in  these  cases  equation  (1)  would  have 
to  be  used,  with  its  non-linear  terms.  The  calculations  cannot  be  carried  out  using  a 
conventional  finite-difference  formulation,  because,  in  essence,  the  system  becomes  non¬ 
isotropic  almost  immediately,  due  to  the  rapidly  varying  diffusivity  coefficient 
In  this  paper,  an  alternative  finite  element  method  for  solving  problems  with  high 
concentration  dependence  in  is  described.  The  problem  requires  a  re-statement  of 
the  domain  properties  at  every  time  step.  The  calculation  of  concentration-dependent 
diffusivity  was  carried  out  using  a  cubic  spline  model.  The  paper  discusses  accuracy  vs. 
computational  effort,  and  also  a  possible  application  in  the  experimental  measurement 
of  diffusivities. 
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F  van  KEULEN 

A  family  of  triangular  plate  bending  elements,  as  a  basis  for  geometrically 
nonlinear  shallow  shell  elements 

A  class  of  shell  elements  can  be  formulated  as  a  combination  of  plate  bending  and 
membrane  elements.  Generally,  this  results  in  facet  elements.  Some  well-known 
examples  are  the  facet  elements  based  on  the  constant  bending  moment  triangle,  the 
Discrete  Kirchhoff  Triangle  and  on  the  Discrete  Shear  Triangle.  In  order  to  account  for 
initially  curved  geometries  and  to  improve  numerical  efficiency,  shallow  shell  terms  can 
be  considered.  This  results  in  shallow  shell  elements,  which  were  originally  facet 
elements. 

Here  our  attention  is  concentrated  mainly  on  a  family  of  plate  bending  elements  which 
is  very  suitabale  as  a  basis  of  the  formulation  of  shallow  shells.  This  family  consists  of 
low-order  as  well  as  higher-order  plate  bending  elements.  The  family  has  members  that 
are  based  on  the  Kirchhoff-Love  assumptions  and  members  that  are  based  on  first  order 
shear  deformation  theory.  The  independent  kinematic  degrees  of  freedom  are  the  nodal 
displacements,  the  rotations  about  the  element  sides  and  additionally,  the  transverse 
shear  strains  along  the  element  sides. 

Some  numerical  examples  are  shown,  which  demonstrate  the  features  of  the  plate 
bending  elements  under  consideration. 

Finally,  the  way  towards  geometrically  nonlinear  shallow  shell  elements  based  on  this 
family  of  plate  bending  elements  is  outlined.  The  formulation  of  efficient  shallow  shell 
elements  is  not  completely  straightforward.  In  particular,  the  higher-order  elements  will 
require  more  attention,  in  order  to  avoid  membrane  locking  and  to  achieve  a  correct 
description  of  rigid  body  motions. 


H  van  LENGEN,  B  MAY,  F-G  BUCHOLZ  and  H  A  RICHARD 

A  substructured  elastic-plastic  fracture  analysis  programme  for 
parallel  processing  on  transputer  networks 

In  recent  years  the  Finite-Element-Method  became  a  powerful  and  excellent 
computational  tool  for  the  analysis  of  complex  problems  in  engineering  science.  In  order 
t  >  enable  the  solution  of  large  problems  on  low  cost  computers  with  small  high  speed 
memory  (RAM)  and  slow  peripheral  storage  media  (Hard  disk)  special  techniques  have 
been  developed.  The  two  most  important  of  these  approaches  are  the  substructuring 
techniques  and  the  frontal  solution  method  / 1,2/.  Using  substructuring  techniques  means 
to  subdivide  the  whole  FE-Model  into  smaller  sub-regions  with  two  classes  of  variables: 
the  internal  variables  and  the  external  variables  related  to  the  nodes  on  main  net  level. 
Using  such  substructured  models  most  steps  of  the  FE-c«lculation  can  be  done  without 
synchronisation  in  parallel  on  different  processors.  The  implicit  parallelism  of  this 
method  in  combination  with  the  frontal  solution  method  provides  an  excellent  basis  for 
the  utilisation  of  modern  hardware  developments  for  massive  parallel  computing  [3]. 

This  paper  presents  a  finite  element  programme  for  complex  elastic-plastic  fracture 
analysis  on  transputer  networks.  In  those  networks  each  transputer  has  access  to  4-16Mb 
RAM-memoiy  and  is  connected  through  a  maximum  four  links  to  other  transputers  of 
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the  network.  The  parallel  version  of  PLASTPARA  has  been  implemented  under  the 
HELIOS  operating  system  both  on  a  small  PC-host  based  network  and  on  the 
Supercluster  SC320  of  PC1  (Paderborn  Centre  of  Parallel  Computing).  The  algorithms 
are  designed  and  optimised  with  respect  to  the  distributed-memory  model  of  transputer 
hardware,  so  that  the  required  I/O  transfer  between  the  processors  and  to  the  Hard  disk 
is  minimised.  All  sources  are  written  in  FORTRAN. 

The  programme  PLASTPARA  is  split  into  a  master  programme  and  several  worker 
programmes  running  on  different  transputers.  So  every  substructure  can  be  calculated 
in  parallel  up  to  the  step  of  reducing  all  variables  of  the  substructures  to  the  external 
ones  by  the  frontal  solution  method.  In  the  next  step  the  master  programme  will 
assemble  and  solve  these  subequations  on  one  transputer  and  compute  the  displacements 
on  main  net  level.  In  the  last  step  the  internal  displacements  are  calculated  by  parallel 
backward  substitution  on  subnet  level.  By  this  method  a  high  degree  of  parallel 
processing  is  already  achieved,  and  will  be  improved  by  a  parallel  main-net  solver  in 
future. 

As  an  example  for  the  efficency  of  this  approach  2D  and  3D  FE-analyses  will  be 
discussed.  The  effect  of  the  number  of  substructures  and  transputers  will  be  given  for 
both  implementations  of  PLASTPARA. 
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M  VANMAELE 

On  an  external  finite  element  method  for  a  second-order  eigenvalue  problem  on 
a  concave  2D-domain  with  Dirichlet  boundary  conditions 

We  consider  a  finite  element  method  for  a  second-order  elliptic  eigenvalue  problem  in 
a  concave  bounded  domain  OCR2  with  sufficiently  regular  boundary  dft  ,  where 
homogeneous  Dirichlet  boundary  conditions  are  imposed.  In  the  method  a  'variational 
crime’  is  committed,  the  approximate  eigenvalue  problem  being  formulated  on  a  domain 
ft*  <t  ft.  Hence  the  finite  element  approximation  space  Vk  <t  V,  V  =  H£(ft). 

We  restrict  ourselves  to  a  linear  triangular  finite  element  mesh  (b  mesh  parameter). 
The  convergence  of  the  method  and  optimal  error  estimates  for  both  the  eigenvalues  and 
the  eigenfunctions  are  obtained  in  [1].  The  proofs  of  the  estimates  in  that  paper  rest 
heavily  upon  the  0(h2)-estimate  for  the  error  of  the  elliptic  projection  in  the  L,(fl*)- 
norm,  see  [1,  Theorem  4.8].  In  the  present  paper  we  prove  this  estimate  in  an 
alternative  way,  proceeding  similarly  as  in  the  proof  of  the  well-known  Aubin-Nitscbe 
trick.  Moreover  in  this  paper  we  obtain  an  optimal,  0(bl), estimate  (from  above)  for  the 
eigenvalues  under  weaker  conditions  on  dft  than  in  [1,  Theorem  5.5], 


66 

1.  M  Vanroae )e  and  A  2,eniSek,  External  finite  element  approximations  of 
eigenvalue  problems,  RAIRO  SfAN  (accepted) 


K  P  WANG  and  J  C  BRUCH,  Jr 

An  efficient  iterative  parallel  finite  element  computational  method 

An  implementation  of  using  parallel  computation  along  with  adaptive  mesb  finite 
element  analysis  will  be  discussed.  In  the  h-refinementof  the  finite  element  analysis,  the 
number  of  nodes  and  elements  will  be  increased  after  a  finite  element  mesh  has  been 
refined.  The  computation  load  for  the  finite  element  system  also  increases.  This 
situation  can  be  resolved  by  using  a  multi-processor  computer  so  that  the  computation 
load  can  be  shared  by  many  processors  after  an  appropriate  load  balancing  operation. 
A  very  efficient  parallel  iterative  scheme  for  the  finite  element  system  was  utilized  to 
meet  the  requirements  necessary  to  solve  the  system  on  a  multi-processor 
(parallel/concurrent)  computer.  In  addition,  different  refinement  schemes  were 
employed  to  investigate  the  effectiveness  of  the  adaptive  mesh  finite  element  analysis 
procedure.  The  model  problem  considered  was  a  free  surface  seepage  problem 
formulated  using  a  fixed  domain  method.  This  imposed  a  restriction  on  the  numerical 
iterative  approach.  Numerical  results  were  obtained  on  an  iPSC/2  Hypercube  concurrent 
computer. 


M  K  WARBY  and  J  R  WHITEMAN 

! 

The  computational  modelling  of  the  thermoforming  process  for  the  creation  of 
axisymmetric  container  structures 

This  paper  is  concerned  with  the  computational  modelling  of  the  thermofonning  process 
which  is  a  process  used  in  the  packaging  industry  to  create  polymeric  container  structures 
by  forcing  thin  sheets  into  moulds.  In  our  model  only  axisymmetric  geometries  are 
considered  and  the  sheet,  although  generally  multilayered,  deforms  to  a  reasonable 
approximation  like  a  membrane.  When  the  containers  to  be  formed  are  'shallow”  it  is 
usually  satisfactory  to  use  pressure  alone  as  the  forcing  action  and  this  case  has  been 
considered  by  several  authors  in  the  recent  literature.  In  this  context,  the  term 
satisfactory  usually  refers  to  the  thickness  distribution  of  the  formed  container  (pot).  In 
our  model  the  process  starts  with  a  sheet  of  uniform  thickness  and  as  the  deformation 
proceeds,  the  thickness  distribution  becomes  non  uniform.  When  the  thickness 
distribution  becomes  highly  non  uniform  or  more  specifically  when  any  part  of  the  sheet 
becomes  too  thin  then  we  would  describe  the  container  as  being  unsatisfactory.  This  is 
the  situation  that  occurs  when  pressure  alone  is  used  and  the  pot  is  deep.  In  the 
manufacturing  process  deep  pots  are  instead  created  by  the  combined  actions  of  a  rigid 
plug  and  an  applied  pressure.  The  use  of  the  plug  leads  to  a  different  deformation 
history  before  contact  is  made  with  the  mould  than  is  the  case  when  pressure  alone  is 
used  and,  since  the  material  is  viscoelastic,  this  has  the  potential  of  leading  to  a 
thickness  distribution  which  is  different  from  that  obtained  with  pressure  alone. 
However  the  main  reason  for  the  more  satisfactory  distribution  of  the  thickness  being 
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obtained  is  as  a  result  of  the  slight  cooling  effect  that  the  plug  has  on  the  sheet 
Specifically  the  plug,  which  is  cooler  than  the  sheet,  has  the  effect  of  slightly  cooling  and 
hence  slightly  stiffening  the  part  of  the  sheet  to  which  it  makes  contact  The  purpose 
of  constructing  a  computational  model  is  to  be  able  to  quantify  this  effect  for  different 
plug  shapes,  depths  and  temperatures.  This  requires  a  model  which  takes  account  of  the 
following  effects. 

The  model  must  take  account  of  the  large  deformation  and  hence  the  equations  are 
nonlinear.  There  is  a  contact  problem  as  the  sheet  slides  on  the  plug  and  for  the  part 
of  the  sheet  in  contact  with  the  plug  there  is  the  nonisothermal  cooling  effect  When  the 
sheet  is  inflated  against  the  mould  there  is  a  simpler  contact  problem  as  the  sheet  sticks 
to  the  mould.  For  the  constitutive  model  for  the  polymers  only  standard  Mooney  Rivlin 
elastic  and  simple  viscoelastic  relations  have  been  considered  as  no  more  better 
experimentally  based  information  was  available.  Fortunately  the  model  involved  only 
one  space  dimension  and  thus  the  finite  element  discretization  is  straightforward.  In  our 
model  cubic  Hermite  elements  using  a  carefully  chosen  moving  mesh  is  used  in  order  to 
capture  the  main  features  at  relatively  low  cost  A  form  of  Newton’s  method  is  used  at 
each  stage  to  solve  the  nonlinear  system  of  equations.  Results  will  be  presented  for  the 
intermediate  membrane  profiles  and  the  final  thickness  distribution  predicted  using  the 
model. 


A  J  WATHEN 

Element-by-element  preconditioning 

Element-by-element  preconditioning  methods  were  introduced  by  T  J  R  Hughes  and 
co-workers  in  1983  as  efficient  and  highfy  vectorisablc  techniques  which  arise  out  of  and 
are  highly  consistent  with  natural  finite  element  datastructures. 

In  this  talk  we  present  a  description  and  analysis  of  such  *EBE'  preconditioning 
methods  for  self-adjoint  elliptic  partial  differential  equations  in  2-  and  3-dimen$ional 
domains.  As  well  as  analysis  for  model  Poisson  problems,  we  present  results  for 
problems  with  discontinuous  variable  coefficients  and  for  constant  coefficient  problems 
discretised  using  elements  with  large  aspect  ratios. 

Some  of  this  work  is  joint  with  Han-Chow  Lee. 


J  WEISZ 

A  posteriori  error  indicator  and  estimator  for  time-space  FE  discretization 
of  parabolic  problems 

In  the  area  of  adaptive  numerical  methods  for  parabolic  partial  differential  equations 
there  are  two  traditional  directions:  In  the  first  one  the  equation  is  discretized  in  space 
and  then  the  corresponding  system  of  ordinary  differential  equations  is  solved  by  usual 
methods  using  the  adaptive  control  of  the  time  step.  The  second  possibility  is  to 
discretize  the  equation  in  time  and  then  to  solve  the  corresponding  elliptic  equations 
adaptively. 
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The  simultaneous  time-space  discretization  [1)  is  a  generalization  of  the  above 
methods.  The  equation  is  discretized  in  time  and  space  simultaneously  using  time-space 
finite  element  spaces.  We  obtain  an  equation  which  is  'elliptic*  and  finite-dimensional. 
Our  aim  is  to  derive  a-posteriori  error  indicators  and  estimators  for  such  discretization 
of  the  linear  parabolic  initial-boundary  value  problem.  Using  the  technique  of  [2]  the 
error  of  the  finite  element  approximation  is  estimated  in  terms  of  the  discretization 
parameter,  data  of  the  problem  and  the  approximate  solution. 

1.  Axelsson,  O,  Maubach,  J,  A  time-space  finite  element  method  for  nonlinear 
convection  diffusion  problems,  NMFM  30,  Vieweg  Braunschweig  1990,  6-23. 

2.  Baranger,  J,  El  Amri,  H,  Estimateurs  a  posteriori  d'erreur  pour  le  calcul  adaptif 
d’ecoulements  quasi-newtoniens,  M2AN  25  (1991),  31-48. 


H  WERN 

Finite-element  solutions  for  mechanical  drilling  methods 

The  principle  for  evaluating  the  unknown  depth  distributions  of  residual  stresses  from 
drilling  methods  relaxed  strain  data  is  viewed  as  an  Integra)  Method  [1-3].  The 
mechanical  hole  drilling  [4,5]  as  well  as  the  ring-core-method  [6,7]  are  analyzed  using 
finite  elements.  The  influence  of  the  notch  sensitivity  is  discussed  for  different  hole-, 
slot-  and  sample  geometries.  Temperature  effects  caused  by  the  drilling  procedure  are 
considered.  From  the  calculated  strain  relaxation  data,  possible  integral  operators  for 
the  relaxation  process  are  discussed. 
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N-E  W1BERG  and  F  ABDULWAHAB 

Improved  patch  stress  recovery  procedure  by  equilibrium  and  boundary  conditions 

The  cornerstone  of  an  adaptive  FE-analysis  is  a  cheap  and  accurate  calculation  of  the 
local  and  global  distribution  of  the  enor,  so  that  successively  improved  solutions  can  be 
obtained.  The  Zienluewicz-Zhu  error  estimate  is  undoubtedly  the  most  practical  a 
posteriori  error  estimator  capable  of  estimating  both  local  and  global  error  of 
discretization.  This  error  estimator,  for  instance  in  the  energy  norm,  is  asymptotically 
exact  if  the  recovered  derivatives  are  superconvergent  Very  recently  two  important 
recovery  procedures  that  recover  superconvergent  gradients  have  been  developed. 
Namely  the  Superconvergent  Patch  Recovery  (developed  by  Zienkiewicz  and  Zhu)  and 
patch  recovery  based  on  superconvergence  and  equilibrium  (developed  by  the  authors). 
Both  methods  are  based  on  local  patches  of  elements.  Their  performance  is  excellent 
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for  recovering  derivatives  for  points  in  inner  patches.  Generally  speaking  the  quality  of 
recovered  derivatives  at  or  near  boundaries  is  however  inferior  compared  to  that  at 
interior  points. 

Both  methods  do  not  take  into  consideration  the  boundary  conditions  which  can  be 
either  prescribed  displacements  or  prescribed  tractions.  In  this  paper,  a  postprocessing 
technique  for  points  on  or  near  boundaries  is  described.  We  introduce  a  weighted  least- 
squares  polynomial  fit  in  an  attempt  to  'force'  the  recovered  derivative  field  to  satisfy  the 
prescribed  boundary  conditions.  In  the  literature,  the  least  square  methods  are  deemed 
to  be  insensitive  to  different  weighting  functions  over  a  wide  range.  In  spite  of  this,  we 
assume  that  higher  values  of  weighting  functions  for  prescribed  boundary  conditions 
would  improve  the  recovered  derivatives  near  boundaries.  Several  numerical  examples 
illustrate  the  very  good  performance  of  the  proposed  recovery  technique  at  or  near 
boundaries. 

J  R  WILLIAMS  and  K  AMARATUNGA 
Matrix  and  image  decomposition  using  wavelets 

This  paper  demonstrates  bow  the  wavelet  transform  may  be  advantageously  applied  to 
matrix  and  image  data.  The  resulting  decomposition  yields  a  reduced  set  of  data  which 
retains  the  essential  properties  of  the  original  data.  As  a  result,  the  wavelet 
decomposition  enables  us  to  capture  the  behaviour  of  a  physical  system  in  a  greatly 
simplified  mathematical  model. 
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Figure  1:  Eigenvalues  of  a  stiffness  matrix  and  wavelet  reduced  matrices  using  (a)  D6 
and  (b)  D20 

Figure  1  compares  the  eigenvalues  of  the  stiffness  matrix  for  a  spring-mass  system  with 
those  of  the  average  matrices  of  ‘/4  and  Vi6  the  size  derived  using  the  Daubechies  D6  and 
D20  wavelet  systems.  The  figure  shows  that  good  estimates  of  the  lower  eigenvalues  of 
the  stiffness  matrix  may  be  obtained  from  the  averaged  matrices.  This  result  has 
significant  implications  in  modal  analysis,  where  the  lower  eigenfrequenciesare  often  the 
most  critical  of  all. 

Similarly,  in  the  wavelet-Galerkin  method  for  solving  PDE’t,  the  resulting  wavelet 
differential  operator  matrix  may  be  decomposed  using  the  wavelet  transform,  leading  to 
an  averaged  differential  operator  matrix  and  a  reduced  set  of  equations.  By  solving  the 
reduced  set  of  equations,  a  good  estimate  of  the  solution  may  be  obtained  with  relatively 
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little  work. 

The  averaging  property  of  wavelets  is  attractive  in  itself  for  the  rapid  solution  of 
engineering  problems.  The  appeal  of  wavelets  is  greatly  enhanced  by  their  hierarchical 
properties,  which  permits  the  subsequent  addition  of  detail  to  the  solution,  thereby 
providing  the  capability  for  progressive  refinement  of  the  solution. 


ZPWU  and  D  PHILLIPS 

A  general  bond-slip  formulation  for  embedded  reinforcing  bar  elements. 

This  is  a  generalized  formulation  developed  for  isoparametric  elements  including 
embedded  reinforcing  bars  and  bond-slip  effects,  in  which  the  representation  of  the 
bond-slip  behaviour  is  evaluated  by  prescribing  bond  characteristics  at  a  set  of  artificial 
"nodes”  introduced  along  the  bar.  The  embedded  bar  element  and  the  parent  concrete 
element  are  now  general  finite  elements  removing  all  the  restraints  from  the 
conventional  approaches.  The  condition  of  compatibility  for  the  system  is  obtained  by 
only  condensing  the  middle  nodes  of  the  bond-slip  while  the  equilibrium  condition  is 
satisfied  by  assembling  the  stiffness  matrix  of  the  reinforcing  bar  and  bond  slip  into  that 
of  the  parent  concrete  element  at  element  level.  Several  examples  are  employed  to 
illustrate  the  performance  of  the  formulation,  i.e.  pull-out  test,  transfer  test  and  beam- 
column  joints.  The  numerical  results  are  compared  with  experimental  ones  where  good 
agreements  have  been  achieved. 
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